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Preface 


The  purpose  of  this  study  was  to  apply  the  boundary 
element  method  (BEM)  to  two  dimensional  fracture  mechanics 
problems,  and  to  use  the  BEM  to  analyze  the  interference  effects 
of  holes  on  cracks  through  a  parametric  study  of  a  two  hole 
tension  strip.  The  study  analyzed  the  effect  of  hole  diameter, 
pitch  and  crack  length.  The  results  of  the  study  were  to  be 
applied  to  a  sample  crack  growth  analysis  to  display  the  use  of 
the  boundary  element  method  in  conventional  aircraft  damage 
tolerance  analysis. 

The  analysis  of  classical  fracture  problems  showed 
excellent  results,  and  the  comparisons  to  different  finite 
element  methods  were  also  very  good. 
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enabling  me  to  complete  this  thesis. 
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Abstract 


This  investigation  analyzes  a  crack  emanating  from 
one  hole,  and  approaching  a  second  hole,  in  a  two  hole 
tension  strip  with  finite  boundaries  using  the  Boundary 
Element  Method.  The  study  included  the  effects  of  varying  the 
hole  diameter,  hole  separation  and  the  length  of  crack.  The 
final  results  were  plotted  as  a  function  of  the  geometric 
correction  factor  0,  which  can  be  presented  as  a  family  of 
curves.  An  example  damage  tolerance  analysis  is  presented 
with  the  3  curves  being  incorporated  into  a  0  look-up  table 
as  used  in  the  NASA/FLAGRO  fatigue  crack  growth  program.  This 
technique  is  acceptable  in  most  fatigue  crack  growth  programs 
now  used  in  the  aircraft  industry  to  ensure  aircraft 
structural  integrity. 

Several  classic  fracture  mechanics  problems  are 
analyzed,  and  computational  efficiency  as  compared  to 
conventional  finite  element  techniques  is  investigated. 
Agreement  with  analytic  solutions  as  well  as  other  numerical 
methods  (finite  element)  is  excellent.  The  computation 
efficiency  was  shown  to  an  improvement  over  existing 


methods . 


Introduction 


I . 


The  present  day  acceptance  of  a  fracture  mechanics 
based  aircraft  damage  tolerance  criteria  is  based  on  the  work 
done  in  the  late  sixties  and  early  seventies,  credited  to  Mr 
Charles  Tiffany  and  Dr  John  Lincoln  [16].  Dr  Lincoln  gives  an 
excellant  review  of  the  Air  Force  Damage  Tolerance  experience 
in  reference  [25].  Within  the  last  decade,  the  Air  Force  has 
placed  increased  emphasis  on  the  fracture  mechanics 
life-cycle  structural  integrity  of  its'  manned  aircraft 
systems  [6].  The  original  implementation  of  an  Aircraft 
Structural  Integrity  Program  (ASIP)  was  in  1959  and  the 
catastrophic  events  leading  up  to  it  are  well  documented  [2]. 
The  most  significant  event  being  the  B-47  fatigue  failures 
which  crippled  the  Strategic  Air  Command  at  a  time  of  extreme 
world  tension.  The  B-47  showed  that  modern  aircraft  could  no 
longer  be  designed  solely  for  static  strength.  This  1959  ASIP 
involved  the  "Safe-Life"  concept  revolving  around  classical 
"Fatigue"  analysis.  To  account  for  the  significant  "scatter 
factor"  associated  with  fatigue  testing,  a  scatter  factor  of 
"four"  was  established  whereby  an  aircraft  designed  for  a 
4000  hour  service  life  must  be  analyzed  and  tested  for  16000 
hours  of  service  life.  The  F-111  was  such  an  aircraft  and  was 
tested  successfully  for  16000  hours.  However,  in  December  of 
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1969  an  F-111  with  approximately  100  flight  hours  crashed  at 
Nellis  AFB,  Utah. 

The  cause  of  this  crash  was  a  manufacturing  defect  in 
the  wing  pivot  fitting  that  was  undetected  by  inspections. 
Also,  the  KC-135  aircraft  was  judged  to  have  a  Safe-Life  of 
13000  hours,  yet  service  experience  had  detected  fourteen 
cases  of  unstable  cracking  in  the  lower  wing  skins  at  flight 
times  ranging  from  1800  to  5000  hours.  An  F-5  which  was 
judged  to  have  a  Safe-Life  of  4000  hours  was  lost  at  Williams 
AFB,  Arizona,  with  approximately  1900  hours.  From  these  and 
other  cases  it  was  apparent  that  the  Safe-Life  methodology 
had  not  precluded  the  use  of  "brittle"  materials  and  "rogue" 
manufacturing  or  service  induced  defects  that  could  lead  to 
premature  failure.  The  re-^ult  was  the  implementation  of  a 
"Fracture  Mechanics"  based  "Damage  Tolerance"  approach  to 
structural  integrity  (5J. 

The  Damage  Tolerance  approach  relies  on  fatigue  crack 
growth  calculations  to  establish  the  time  interval  required 
to  grow  a  crack  from  an  initial  size  (usually  the  maximum 
flaw  undetectable  with  current  NDI  techniques)  to  the 
critical  crack  length  which  denotes  the  onset  ot  unstable 
crack  growth.  The  crack  growth  equations  are  a  function  of 
the  local  change  in  Stress  Intensity,  K,  as  a  stress  cycle  is 
applied.  All  of  the  current  fatigue  crack  growth  codes  in  use 
by  industry  have  "canned"  subroutines  to  calculate  K  for 
classical  crack  configurations.  However,  the  practicing 
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engineer  is  frequently  faced  with  design  details  which  are 
not  represented  by  the  conventional  solutions.  In  this  case 
the  solution  for  K  must  be  analyzed  independently  of  the 
crack  growth  code,  with  the  results  placed  in  a  look-up  table 
as  a  function  of  crack  length. 

The  most  common  method  of  independent  analysis  of 
Stress  Intensity  Factors  is  the  finite  element  method.  This 
numerical  technique  is  extremely  flexible  in  its'  ability  to 
analyze  a  wide  range  of  problems.  However,  finite  element 
models  require  the  descretization  of  the  body  being  studied, 
with  a  gradual  refinement  towards  the  crack  tip.  This  is  very 
expensive  in  both  computer  time  and  man  hours.  Alternative 
solution  techniques  are  always  being  sought  to  increase  the 
efficiency  of  these  Stress  Intensity  Factor  analyses,  and 
this  thesis  will  investigate  the  possible  advantages  of  using 
the  Boundary  Element  Technique. 

The  Boundary  Element  Method  is  based  on  a  singular 
solution  which  represents  the  analysis  of  a  segment  of  the 
boundary  of  the  body  being  analyzed.  The  values  of  the 
boundary  conditions  are  known,  and  the  solution  calculates 
the  results  for  the  rest  of  the  body.  The  singular  solution 
will  satisfy  the  governing  differential  equations  exactly, 
and  the  user  will  approximate  the  boundary  conditions. 

Some  of  the  earliest  uses  of  BEM  were  in  1963  by 
Jaswon  and  Ponter  [10],  Jaswon  [9],  and  Symm  [22]  concerning 
potential  problems.  The  first  elasticity  application  was  in 
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1967  by  Rizzo  [17).  Since  then  the  method  has  been  used  in  a 
wide  variety  of  applications  as  documented  by  Mackerle  and 
Andersson  [12].  Examples  of  the  BEM  applied  to  fracture 
mechanics  can  be  found  in  papers  by  Snyder  and  Cruse  [20]  and 
Rizzo  and  Shippy  [18]. 

The  purpose  of  this  research  is  to  investigate  a 
crack  emanating  from  a  hole  towards  another  hole  in  a  two 
hole  tension  strip  with  finite  boundaries.  The  basic  BEM 
technique  will  be  verified  on  similar  problems,  either 
classical  or  by  finite  element  methods.  A  parametric  case 
study  of  the  two  hole  tension  strip  was  conducted  varying  the 
hole  diameter,  hole  separation  and  crack  length  to  create  a 
family  of  Stress  Intensity  Factor  data  curves  suitable  for 
fatigue  crack  growth  analysis.  A  sample  damage  tolerance 
analysis  using  the  results  of  the  parametric  stress  intensity 
study  is  shown. 

After  a  theoretical  development  of  the  Fictitious 
Stress  BEM  technique,  comparison  solutions  to  a  few  classic 
fracture  problems  are  presented,  as  well  as  a  comparison  to 
finite  element  solutions.  The  solutions  to  the  two  hole 
tension  strip  parametric  study  are  presented,  followed  by  the 
example  damage  tolerance  analysis. 
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II .  Theoretical  Discussion 

The  analytical  method  used  for  this  study  was  the 
Boundary  Element  Method.  This  method  relies  on  a  singular 
solution  representing  the  analysis  of  a  segment  of  the 
boundary  of  the  body  being  analyzed.  The  boundary  conditions 
are  known,  and  the  solution  analyzes  the  impact  on  the 
remainder  of  the  body.  The  complete  solution  requires  all  of 
the  boundary  segments  be  solved  simultaneously  and  include 
the  effects  of  the  boundary  segments  on  each  other.  The 
technique  used  in  the  analysis  performed  here  is  the 
"Fictitious  Stress"  method  as  outlined  by  Crouch  and 
Starfield  [4]. 

A.  Fictitious  Stress  Method 

The  Fictitious  Stress  method  utilizes  the  plane 
strain  version  of  Kelvin's  problem  [21]  as  the  basic  singular 
solution.  Kelvin's  problem  is  a  point  load  in  an  infinite 
domain  while  the  plane  strain  version  is  a  line  of 
concentrated  force. 

As  shown  in  Figure  1,  the  plane  strain  version  of 

Kelvin's  problem  will  involve  a  line  of  force  F  along  the  z 

direction.  The  components  F  and  F  have  units  of  force  per 

X  y 

unit  depth.  Kelvin  showed  that  a  harmonic  function  g(x,y)  was 

2  2 

a  solution  to  the  biharmonic  equation  (7  7  -0)  such  that 
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g(x,y)  -  [-ln(x^  +  ) ^^^ ]/( 4n( 1-v) ] 


(1) 


where  v  is  Poisson’s  ratio. 

Kelvins  solution  for  displacements  components,  and 
Uy,  are  expressed  as 

Ux  -  (  Fj^/2G)  [  (  3-4v)g  -  +  (  Fy/2G )  [ -yg  )  (2) 

"  (F„/2G)[  (3-4v)g  -  yg  ]  +  (F^/2G){-xg 
y  y  ^  y  ^  » y 

where  F  and  F  are  the  components  of  the  applied  point  load 
X  y 

F,  and  G  is  the  material  shear  modulus,  and 


g  -  g(x,y)  (3) 

g  ,  -  [-l/4ii(l-v)]tx/(x2  +  y^)J 

g^y  "  [-l/4ll(l-v)  ][y/(x2  +  y2)) 

g  -  (l/4ii(l-v)  I[2xy/(x2  +  y2)2j 

f  *y 

9  XX  -  -5  vv  -  [1/4r(1-v)]1(x^  -  y^)/(x2  +  y^)^] 
f  ^  y  y 

where  g  ^  denotes  the  partial  differentiation  of  the  function 
g(x,y)  with  respect  to  y  (3g/3y)  and  g  ^  denotes  the  partial 
differentiation  of  g  with  respect  to  x  (3g/3x).  Using  the 
stress-strain  relations  produces  the  stress  results  as 
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(4) 


®xx 

*9,xx' 

+ 

Fyl2M9_y  -  y9,^^l 

"yy  “ 

y9,i,yi 

+ 

-  ''9,yy> 

"xy  “ 

F,Ul-2v)9_y  - 

+ 

Fyl(l-2v)g  -  yg 

.  xy ' 


where  <t  <r  „  are  the  components  of  stress  in  the  x  and  y 
XX  yy 

directions  and  a  is  the  shear  stress.  The  stresses  in  Eq 

xy 

(4)  satisfy  the  equations  of  equilibrium,  and  are  singular  at 
the  point  x*y«0.  Timoshenko  and  Goodier  showed  that  these 
stresses  correspond  to  a  line  of  concentrated  force  at  the 
origin  [ 24 ] . 

To  facilitate  the  transition  of  Kelvin's  solution  into  a 

form  usable  in  a  numerical  technique,  we  consider  the  problem 

of  tractions  t  -  P  and  t  -  P  applied  to  a  line  segment 

XX  y  y 

an  infinite  elastic  solid.  Kelvin's  solution  is 
integrated  over  a  line  segment  of  length  2a  as  shown  in 
Figure  2.  If  we  consider  a  small  segment  of  the  line,  de,  the 
force  F  becomes 


Fi(e)  -  P^de  i-x,y  (5) 

We  will  assume  a  constant  traction  thus  the  new  harmonic 
function  to  satisfy  the  biharmonic  equation,  f(x,y),  as  shown 
in  [4],  can  be  expressed  in  terms  of  g(x,y)  as 
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Figure 


f(x,y)  =/■  g(x  -  8,y)  cfe 


2.  Integration  of  Kelvin's  Problem 
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f(x,y) 


(6) 


g^J®g(x-e,y)dc 

-  [-l/4n( 1-v) 1 (y{arctan(y/(x-a) ) 

-  arctan(y/x+a) } 

+  (x+a)ln{(x+a)^  + 

-  (x-a)ln{(x-a)^  + 

Following  the  procedure  outlined  previously  in  the 
presentation  of  Kelvin's  problem,  the  displacements  due  to 
the  line  of  concentrated  force  per  unit  depth,  Fj^{c),  are 

Ux  -  (Px/2G)[(3-4v)f  +  yf  y)  +  ( Py/2G) ( -yf  ^ 1  (7) 

Uy  -  (Py/2G)((3-4v)f  +  yf^y)  +  < Px/2G) t -yf ^ x ^ 

and  the  stresses  become 

%x  -  *  y'.xyl  *  >f‘,yy>  <®> 

’yy  -  ’■x>-l'l-2''>',x  -  y'.xy'  *  Pyl 2( l-^)* ,y  +  yf,yyl 
'xy  -  *  >f',yy>  *  'y' <l-2''>*,x  *  y‘,xy> 

The  derivatives  of  f  are  given  as 

f,x  -  l/[4n(l-v)l(ln{(x-a)2  +  (9) 

-  ln{(x+a)^  + 

f  y  ■  -l/[  4»i(  l-v)  1  [  arctan{y/(  x-a ) }  -  arctan{y/(  x+a ) }  ] 

f  -  l/[4ii(l-v)  ][y/{(x-a)2  +  y^)  -  y/{(x+a)^  +  y^}] 

f 
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^,xx  ”  ^,yy 

-  l/(4ii(l-v)  ][  (x-a)/{(x-a)^  +  y^} 

-  (x+a)/{(x+a)^  +  y^}J 

It  is  important  to  note  that  the  stress  solutions  are  not 
defined  for  x-+a,  or  y-0.  To  investigate  this  further,  it  is 
necessary  to  consider  the  the  stress  tensor  along  the  line 
y«0.  Evaluating  Eqs  (8)  and  (9)  for  y-0  yields  the  stresses 

-  -(3-2v)/[8n(l-v)  ]P^lnt  (x+a)/(x-a)  (10) 

-  2v/( 4n( 1-v) ]Pylimy^Q^(arctan{y/(x-a) } 

-  arctan{y/(x+a ) } ] 

Oyy  -  (1-2v)/{8ji(1-v)  }Pj^ln[  (x+a)/(x-a) 
l/(  2ii)Pylimy^Q^(  arctan{y/(x-a) } 

-  arctan{y/( x+a ) } J 

ffjjy  ”  -l/(2ii)Pj^limy_^Q^(arctan{y/(x-a) } 

-  arctan{y/( x+a ) } 1 

-  (l-2v)/(8n(l-v) lPyln((x+a)/(x-a) 1^ 

where  the  limits  on  y  are  necessary  as  the  arctan  function  is 
multivalued.  The  arctan  functions  in  Eq  (10)  are  interpreted 
to  represent  the  angles,  and  02,  from  the  ends  of  the  line 
segment,  to  an  arbitrary  point  (x,y),  as  shown  in  Figure  3. 
The  values  of  0^^  and  02  are  seen  to  be 

0j^  -  arctan(y/(x-a)  ]  (11) 
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Cl 


©1  =  arctan(y/(x-a)) 
©2  =  arctanV/(x+a)) 


re  3.  Boundary  Element  Geometry 
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arctan[ y/( x+a )  ] 


When  y-0,  it  can  be  seen  that  9  can  be  -n,  0,  or  +n.  By 
examining  the  limit  expression  in  Eq  (10),  we  see  that  the 
three  possible  solutions  are 

limy->o[arctan{y/(x-a) )  -  arctan{y/(x+a ) }  ]  (12) 

-  0  |x|>a,  y-0^  or  y-0_ 

-  +n  |xl<a,  y-0_^ 

-  -n  |x|<a,  y-0_ 

Examining  the  last  two  values  of  Eq  (12),  it  can  be  seen 
the  stress  tensor  is  discontinuous  across  the  line  segment  at 
y«0.  It  is  instructive  to  examine  the  magnitude  of  the 
difference  in  the  stress  tensor  across  y-0  for  |x|<a.  The 
change  in  is 

axx(x,0_)  -  -(3-2v)/l8!i(l-v)]P^ln[(x+a)/(x-a) 

+  PyV/[2(l-v)l 

®xx^*'®+^  -  -(  3-2v)/[8ii(l-v)  iPj^lnl  (x+a)/(x-a) 

-  PyV/I2(l-v)] 

<Txx^*'®-)  -  '^xx^^'®+^  “  PyV/(l-v)  (13b) 


(13) 


(13a) 
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the  change  in  is 


a  (x,0  )  -  (l-2v)/(8ii(l-v)  JP  ln{  (x+a)/(x-a)  ] 

yy  -  A 

+  P/2 


(14) 


a  -  (l-2v)/(8n{l-v) JP  ln( (x+a)/{x-a) 1 

yy  +  A 

-  v 


(14a) 


OyyO^.O-l  -  ^ 


yy 


(14b) 


and  the  change  in  a  is 

xy 


<j^^(x,0_)  -  -(l-2v)/(8rt(l-v)  IP  Inf  (x+a)/(x-a)  J 
xy  "  y 

+  P^/2 


(15) 


(x,0  )  -  -(l-2\>)/[8n(l-v)  ]P  ln[  (x+a)/(x-a)  1' 
xy  +  y 

-  Px/2 


(15a) 


xy 


(15b) 


It  can  be  seen  that  Eqs  (14b)  and  (15b)  indicate  that  the 

stresses  P  and  P  are  the  constant  discontinuities  in  a  and 
X  y 

(Tyy  respectively.  Crouch  and  Starfield  (4)  showed  that  the 
physical  significance  of  the  stresses  P  and  P  could  be 
interpreted  as  imagining  the  line  segment  |x|<0,  y«0  as  a 
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crack  in  an  infinite  elastic  solid.  As  shown  in  Figure  4,  the 
outward  normal  to  the  positive  side  of  the  crack  y*0^  has 
components  n^-(0,-l),  and  the  outward  normal  to  the  negative 
side  y“0_  has  components  n^«(0,l).  With  the  tractions  t^ 
defined  as, 


- 


the  tractions  on  the  two  surfaces  become, 


(16) 


-  -»yy(X,0^) 

t^(X,0.|  -  «^y<X,0.) 
ty<X.O-)  -  'yy**.*!.) 


The  resultant  stresses  obtained  by  adding  the  traction 
components  t^  on  both  sides  of  the  crack.  Substituting  the 
values  of  a  from  Eq  (14b)  and  a  from  Eq  (15b)  into  Eq 

yy  *y 

(17)  yields. 


t^{x,0)  -  (18) 

ty(X,0)  -  Py 


Thus,  the  stresses 
tractions  across  the 


and  Py  represent  the  constant  resultant 
line  segment  |x|<a,  y«0. 


L5 


B,  Numerical  Algorithm 


To  apply  the  fictitious  stress  method  to  a  general 
problem,  it  is  instructive  to  consider  the  case  of  a  hole 
with  a  boundary  C  in  an  infinite  plate.  We  will  let  the  hole 
be  loaded  by  an  outward  pressure  (p)  load  as  depicted  in 
Figure  5. (a).  As  the  hole  boundary  is  otherwise  traction 
free,  it  is  assumed  the  shear  stress  is  zero,  therefore,  the 
known  boundary  conditions  relative  to  the  normal  (n) 
tangential  or  shear  (s)  directions  are 

«T„  -  -p  (19) 

<j  »  0 

s 

We  now  create  a  system  of  N  line  segments,  joined  end  to 
end,  along  a  boundary  C'  as  depicted  in  Figure  5.(b)  that 
represents  the  boundary  C  of  the  hole  in  Figure  5. (a).  Each 
line  segment  i  is  individually  formulated  with  the  fictitious 
stress  solution  for  a  stress  Bach  line  segment  in  this 

example  will  be  of  a  uniform  length  2a.  If  the  length  2a  is 
small  enough  the  boundary  C  will  be  quite  closely  modeled  by 
C'.  The  local  coordinates  n  and  s  as  depicted  in  Figure  5. (a) 
are  relative  to  C  so  they  will  change  depending  on  the 
location  of  the  point  desired.  The  local  coordinates  n  and  s 
in  Figure  5.(b)  will  be  relative  to  each  line  segment  i.  It 
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c 


BEM  Line  Segments 

Figure  5.  Numerical  Method 
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will  then  be  important  to  order  the  line  segments  such  that 
the  local  coordinates  for  C  and  C'  correspond. 

It  is  now  important  to  remember  that  the  line  segments 
are  based  on  Kelvin's  solution  for  a  point  load  in  an 
infinite  solid.  So  each  line  segment,  or  boundary  element,  is 
in  fact  a  line  of  constant  local  stresses  in  an  infinite 
elastic  body,  which  happen  to  coincide  with  the  boundary  C. 
Each  element  will  have  its'  own  applied  stress,  but  each 

element  will  be  affected  by  all  of  the  other  elements.  Using 
the  theory  of  superposition,  if  we  were  to  calculate  the 
stress  at  a  point  in  the  body,  we  would  have  to  sum  all  of 
the  solutions  for  that  point  due  to  all  of  the  elements  i, 
each  element  with  an  applied  stress  P^.  So  to  calculate  the 
final  stresses  ^<r  and  ^<t  at  each  element  i,  at  its' 

S  11 

midpoint,  requires  a  summation  of  the  form 


^  ft  m  3o  1  X  ji 


^  '  j-ll 


]  i-1  to  N 


where  _ ,  and  ^^A _  are  the  boundary 

ss  sn  ns  nn  ■* 

coefficients.  As  an  example,  gives  the  actual  normal 

stress  at  the  midpoint  of  element  i  due  to  the 

application  of  a  unit  normal  stress  to  element  j  (^P^-1). 

Since  the  values  of  and  are  known  from  the 

s  n 

boundary  conditions  of  the  original  problem,  it  remains  for 
the  system  of  applied  stresses  ^P^  and  ^P^^  to  be  solved 
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for  by  assembling  a  system  of  2N  simultaneously  linear 
algebraic  equations  in  as  many  unknowns. 


0 

-P 


ss 


j-1' 


nN 


j-1' 


^nn^^n' 


i-1  to  N 
i-1  to  N 


(21) 


As  described  in  [4],  it  is  important  to  realize  that  the 
stresses  and  are  fictitious,  and  do  not  really  exist. 

They  are  merely  the  system  of  stresses  applied  to  each 
individual  element  along  C'  such  that  the  simultaneous  system 
of  integrated  Kelvin's  solutions  result  in  the  calculation  of 
the  actual  boundary  conditions  of  the  problem  being  analyzed. 


C.  Co-ordinate  Transformation 


The  equations  for  the  transformation  of  each  individual 
elements  local  influence  coefficients  into  a  common  "global” 
system  is  described  in  detail  by  Crouch  and  Starfield  [4].  By 
examining  Figure  6,  we  will  label  the  local  co-ordinate 
system  for  an  arbitrary  element  as  x'  and  y' .  The  element  is 
defined  as  |x'|<a,  y'-O.  The  stresses  applied  to  this  element 
are  P^,  and  P^, . 

The  local  co-ordinate  system  is  produced  with  a 

translation  of  c  in  the  global  x  direction,  c  in  the  global 

X  y 

y  direction,  and  a  rotation  X  about  the  global  z  axis 
(positive  direction  being  counterclockwise).  The  co-ordinate 
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Figure  6.  Line  Segment  of  Arbitrary  Orientation 
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transformation  is  then  as  follows: 


X'  -  (x  -  c  )cosX  +  (y  -  c  )sinX  (22) 

X  y 

y'  -  -(x  -  c  )sinX  +  (y  -  c  )cosX 
X  y 

Substituting  (22)  and  (9)  into  (7)  and  (8)  produces: 

Ux,  -  Pjj,/(2G)I  (3-4v)F^  +  y'Fj]  +  Py,/(  2G)  [-y' F2  ]  (23) 

Uy,  -  Py,/(2G)l  (3-4v)Fj^  -  y'F3]  +  Px,/(  2G)  [-y' F2  ] 


and 


"x'x'  •  Px'^  y'^4^  ^y'l^vFj  -  y'Fg]  (24) 

<^y,y,  -  Pjj,[-(1-2v)P2  -  y'F^l  +  Py,(2(l-v)F3  +  y'Pg] 

«Tx,y,  -  Px,(2(1-v)F3  -  y'Fg]  +  Py,((l-2v)F2  -  y'F^] 

where  the  functions  Fj^...Fg  are  defined  as: 

Fj^  -  f  (x'  ,y' )  (25) 

-  -l/(  4n(  1-v)  Jly'  (arctan(yVlx'-a) ) 

-  arctan(y V(x'+al ) )  -  (x'-a)ln{ { (x'-a)^  + 

+  (x'+a)ln{( (x^+a)^  +  ) 

^2  "  ^  r  X ' 

-  l/[4n(l-v) ] {ln[ (x'-a)^+y'^)^/^  -  ln[ ( x ' +a ) ^+y ' ^ 1 
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-  -1/t  4ii(  1-v)  ]  {arctan(yV(x'-a)  1  -  arctan[yV(x'+a)  ] } 
^4  "  ^,x'y' 

-  l/[4n(l-v) ]{yVt {x'-a)^+y'^]  -  y '/[ (x'+a)^+y'^l ) 

-  1/i 4n(l-v)]{{x'-a)/[(x'-a)^+y'^] 

-  (x'+a)/[(x'+a)^+y'^l} 

To  calculate  the  displacements  and  stresses  at  a 
particular  element  midpoint,  it  is  necessary  to  calculate  x' 
and  y'  as  coordinates  relative  to  the  local  element  location 
and  orientation.  The  calculated  displacements  and  stresses 
from  Eqs  (22)  and  (23)  are  also  in  the  local  x'y'  system. 
Since  this  is  not  convenient,  one  more  transformation  to  Eqs 
(22)  and  (23)  to  compute  the  resulting  displacements  and 
stresses  in  the  global  xy  coordinate  system.  The  relations 
between  the  xy  global  system  and  the  x'y'  local  system  are 

“  Uj^'cosX  -  Uy'sinX  (26) 

Uy  ■  Uy'cosX  -  Uy'sinX 


2 

<^„,.-rCOs  X  -  2a„ sinXcosX  + 

•  2. 
a  ,  , Sin  A 

XX 

X  '  X  X  '  y ' 

y'y' 

2 

®..,..,cos  X  -  20„  sinXcosX  + 

,  .sin^X 

yy 

y 'y '  x'y ' 

X  '  X  ' 
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‘^xy  “  "  Oy,y,  )sinXcosX  +  ff^,y,(cos  X-sin  X) 

Substituting  Eqs  (26)  and  (27)  into  Eqs  (22)  and  (23)  yields 

,/( 2G) [ ( 3-4v) FlcosX  +  y ' ( F2sinX+F3cosX) ]  (28 

+  Py,/(  2G)  (-(  3-4'\))FlsinX  -  y '  (  F2cosX-F3sinX)  1 
u  -  P  ,/(2G) [ ( 3-4v)FlsinX  -  y ' ( F2cosX-F3sinX) ] 

y  ^ 

+  Py,/( 2G) l-( 3-4v)F1cobX  -  y ' ( F2sinX+F3cosX) 1 

■  P^, tF,+2(l-v) (F2cos2X-F3sin2X)  (29 

+  y ' ( F4cos2X+F5sin2X)  I 
+  Py,  tF2-(l-2\))  (F2sin2X+F3cos2X) 

+  y ' ( F4sin2X-F5cos2X) ] 

<t  -  P  ,  lF,-2(l-»)(F2cos2X-F3sin2X) 

-  y' ( F4cos2X+F5sin2X) J 

+  Py, [ F2+( l-2v) ( F2sin2X+F3cos2X) 

-  y ' ( F4sin2X-F5cos2X) ] 

®xy  ”  Pjj,  I  2(  1-v)  (  F2sin2X+F3cos2X) 

-f  y '  (  F4sin2X-F5cos2X)  J 
+  Py , [ ( l-2v) ( F2cos2X-F3sin2X) 

-  y' ( F4cos2X+F5sin2X) J 

As  can  be  seen,  Eqs  (28)  and  (29)  facilitate  the 


computation  of  influence  coefficients  to  express 
displacements  and  stresses  in  terms  of  P  ,  and  P  , . 


D.  Influence  Coefficients 


To  calculate  the  final  influence  coefficients,  one  final 
transformation  of  Eqs  (28)  and  (29)  are  necessary.  Eqs  (28) 
and  (29)  calculate  the  stresses  and  displacements  at  the  i'th 
element  in  the  global  coordinate  system.  We  are  interested  in 
displacements  and  stresses  at  the  midpoint  of  the  i'th 
element  in  i'th  elements  local  coordinate  system,  x',  y',  as 
shown  in  Figure  7.  The  final  transform  is 

X'  -  x'cosY  +  y'sinr  (30) 

y'  -  -x'siny  +  y'cosr 

where  Therefore, 

^u-,  -  ^Ujj,cosY  +  ^Uy,sinY  (31) 

^u-,  -  -^u^,sinr  +  ^Uy,cosY 

and 

^®x'x'“^*^x'x'^°®^‘'^  *  2^ajj,y,sinYCosY  +  ^<yyfy»sin^Y  (32) 

“  2^o^,y,sinYCosY  +  ^ffy,y,cos^Y 
^®x'y'”“^  ^®x'x'"^®y'y' *  ^a^,y,  (cos^Y-sin^Y) 
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By  realizing  that 


j 

j 


n 


i 


i 


u 


n 


n 


(33) 


(34) 


(35) 


and  substituting  Eqs  (31)  and  (32)  into  Eqs  (23)  and  (24) 
produces 

^u  -  /( 2G) ( ( 3-4v)Flcosr  +  y' ( F26inY-F3cosY) ]  (36) 

D  5 

+  ^P^/(  2G)  ( (  3-4v)FlsinY  -  y' ( F2cosY+P’3sinY)  1 
^u  -  ^P  /( 2G) [ -( 3-4^) FlsinY  -  y ' ( F2cosy+F3cosy) J 

11  o 

+  ^P^/( 2G) [ ( 3-4v) FIcosy  +  y' ( F2sinY-F3coSY) 1 

and 

-  ^P„l-2{l-v)(F2sin2Y  -  F3cos2y)  (37) 

s  s 

-  y ' ( F4sin2Y+F5cos2Y) ) 

+  ^P^[ ( l-2v) ( F2cos2Y+F3sin2Y) 

-  y ' ( F4cos2Y-F5sin2Y) 1 
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i 


n 


j 

+ 

+ 


Pgl F2-2 ( 1-v) ( F2cos2y  +  F3sin2Y) 
y' ( F4cos2 Y-F5sin2Y) ] 

^PnC  F^-( l-2v) ( F2sin2Y-F3cos2Y) 
y' ( F4sin2Y+F5cos2Y) 1 


Thus  Eqs  (35)  and  (36)  can  be  expressed  as 


u. 


.N  i  j, 

j-1  ^SS  "s 


j-1  ”sn  '^n 


u 


n 


.  i^^B  ^P  +  .  i^^B  ^P 

j“l  ns  s  j*l  nn  n 


and 
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i 

i 


E*^.  ^P  +  E*^.  ^P 

j»l  ss  s  sn  n 

E^.  /^A  ^P  +  e”.  i^^A  ^P 
j“l  ns  s  nn  n 


(39 


where  etc.,  are  the  final  influence 

ss  ss 

coefficients . 

The  final  matrix  includes  two  sets  of  2N  equations  in  2N 
variables,  one  for  displacements,  one  for  stresses.  However, 
both  sets  of  equations  have  the  fictitious  stresses  ^P^  and 
^P^  as  the  un)inowns.  Therefore,  to  create  a  solvable  system 
of  2n  equations,  of  the  four  boundary  conditions  for  an 
element  i  ( ,  ^u  ,  ^cr  ,  and  ^v_),  only  two  need  be  Icnown 
(one  shear,  one  normal).  The  final  matrix  of  influence 
coefficients  ( 2N  by  2N)  will  consist  of  A's  and  B's  as 
determined  by  the  type  of  boundary  condition  given  for  each 
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element.  Once  the  quantities  and  are  known,  Eqs  (37) 

and  (38)  can  be  used  to  calculate  the  remaining  unknown 
boundary  conditions,  and  influence  coefficients  can  be 
calculated  to  analyze  displacements  and  stresses  at  any  other 
point  in  the  body. 

It  should  be  noted  that  the  resulting  2N  by  2N  matrix  of 
influence  coefficients  is  fully  populated,  as  every  element 
effects  all  other  elements  as  well  as  itself  [4].  This  is  in 
contrast  to  the  banded  stiffness  matrix  produced  by  the 
finite  element  method.  Though  it  will  be  shown  in  this  study 
that  the  boundary  element  method  can  analyze  certain  types  of 
problems  in  far  less  degrees  of  freedom  than  the  finite 
element  method,  the  boundary  element  method  cannot  take 
advantage  of  a  banded  matrix  so  much  of  the  computational 
advantages  are  lost. 

Another  point  of  interest  is  each  elements  "self 
effects".  By  examining  Eq  (6)  we  see  that  the  value  of  the 
integrated  Kelvin's  solution  decreases  with  increasing 
distance  from  the  midpoint  of  an  element.  Therefore,  the 
maximum  value  for  an  influence  coefficient  must  be  for  an 
elements  influence  on  itself.  Crouch  and  Starfield  [4]  show 
that  the  values  of  all  elements  self  effects  are 

'\s  ”  “^n  =  ±0-5  (40) 

“Bgg  -  “  -(3-4v)/[4iiG(1-v)  ](^a)ln(^a)  (41) 
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As  can  be  seen  from  Eqs  (14)  and  (15),  the  stresses  are 
discontinuous  across  an  element.  The  convention  established 
by  Crouch  and  Starfield  dictates  that  "the  boundary  of  a 
finite  body  is  transversed  in  the  clockwise  sense,  whereas 
the  boundary  of  a  cavity  is  traversed  in  the  counterclockwise 

sense".  This  allows  and  ^^A _  to  be  equal  to  0.5 

always . 

E .  Modeling  Considerations 

As  has  historically  been  the  case  with  the  finite  element 
method,  an  engineers  ability  to  "model"  a  problem  correctly 
plays  as  much  a  role  in  the  value  of  the  final  results,  as 
does  the  accuracy  of  the  method  being  used.  The  boundary 
element  method  also  shares  this  characteristic.  Of  particular 
interest  is  the  fact  that  the  user  should  not  calculate 
displacements  or  stresses  for  a  point  "too  close"  to  an 
elements  midpoint  [4].  The  reason  is  that  it  has  been  found 
empirically  that  the  numerical  solution  is  generally 
unreliable  at  points  within  a  circle  of  radius  equal  to  one 
element  length  (2a)  centered  at  the  midpoint  of  a  boundary 
element,  except  at  the  midpoint  itself.  Therefore,  to  obtain 
data  close  to  a  boundary,  the  user  is  forced  to  refine  the 
lengths  of  the  boundary  elements  in  a  gradual  fashion  as  the 
area  of  interest  is  approached.  F.R.  Harris  [8]  developed  a 


I 


3 


LI 


that  can  be 


modeling  technique  for  a  crack  of  length  "a" 
incorporated  into  the  work  done  herein.  His  method  can  be 
stated  as  follows:  The  crack  length  is  divided  into  .50a, 

.25a,  .125a  and  .125a  segments.  The  first  segment,  or  the 
.50a  length  segment  is  divided  into  three  equal  length 
boundary  elements  (element  length  -  .1667a).  The  second 
segment,  or  the  .25a  length  segment,  is  divided  into  three 
equal  length  boundary  elements,  (element  length  -  .0833a). 

The  third  segment  is  divided  into  three  equal  length  boundary 
elements  (element  length  «  .04167a).  The  last  segment  is 
divided  into  25  equal  length  boundary  elements  (element 
length  -  .005a).  By  using  this  method  of  gradual  refinement, 
stresses  can  be  computed  with  reasonable  accuracy  near  the 
area  of  the  singularity  at  the  crack  tip. 

To  model  the  problems  in  this  study,  each  body  was 

modeled  with  a  line  of  symmetry  along  the  line  of  the  crack. 

The  entire  line  of  symmetry  was  modeled  with  boundary 

elements.  The  crack  itself  is  modeled  with  the  F.R.  Harris 

method  of  refinement  [8]  outlined  above.  The  refinement 

scheme  is  mirrored  at  the  crack  tip  both  along  the  crack 

itself,  and  along  the  uncracked  material  directly  in  the  path 

of  the  crack.  The  elements  along  the  non-cracked  boundary 

utilized  enforced  displacement  conditions  (u^»0)  normal  to 

the  line  of  the  crack  and  stress  conditions  tangential  to  the 

crack  (o  -0).  The  crack  surface  itself  was  modeled  as  being 

stress  free  (a^-a^-O). 

n  s 
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Ill .  Boundary  Element  vs  Finite  Element 

An  initial  configuration  of  a  two  hole  tension  strip 
was  analyzed  with  both  the  boundary  element  method  described 
in  this  study,  and  with  the  MSC/NASTRAN  finite  element  code. 
Both  methods  were  used  in  order  to  compare  the  relative 
strengths  and  weaknesses.  Both  method  only  modeled  the  top 
half  of  the  tension  strip  using  symmetry  conditions  as 
enforced  through  restrained  vertical  displacement  along  the 
line  of  symmetry. 

The  geometry  and  material  properties  of  the  problem 
are  illustrated  in  Figure  8.  The  tension  strip  has  two  holes 
of  0.25  inch  diameter.  The  holes  are  separated  by  one  inch. 
The  edge  distances  to  the  holes  are  three  diameters  for  all 
four  sides.  The  initial  crack  length  is  0.1  inch,  and  it  is 
emerging  from  one  hole  and  oriented  towards  the  second  hole. 
The  far  field  tension  stress  is  46  KSI. 

A.  Finite  Element  Method 


There  were  three  MSC/NASTRAN  models  constructed.  This 
was  to  provide  convergence  data.  The  baseline  model  consisted 
of  13,858  degrees  of  freedom.  The  other  two  models  had 
respectively  8,610  and  15842  degrees  of  freedom.  Needless  to 
say  all  of  the  models  were  constructed  with  a  graphic 
pre-processor/model  generator  ( PDA/PATRAN ) .  The  baseline 
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Figure  8.  Tension  Strip  Problem 


33 


model  is  shown  in  Figure  9.  The  course  model  (8,610  DOF)  and 
the  fine  model  (15,842  DOF)  look  identical  to  the  baseline 
model  only  differing  in  the  density  of  the  mesh  in  the 
immediate  vicinity  of  the  crack  tip. 

All  three  models  consisted  of  a  mesh  of  eight  noded 
quadratic  isoparametric  quad  elements  in  the  crack  area.  The 
eight  noded  quad  ( CQUAD8 )  mesh  then  transitions  into  a  four 
noded  quad  (CQUAD4)  element  mesh  to  complete  the  model.  A 
handful  of  six  noded  triangles  (CTRIA6)  were  required  in  the 
transition  region. 

The  entire  models  were  declared  "Surfaces”  as  described 
in  the  MSC/NASTRAN  Users  Manual  [14]  and  interpolated  stresses 
were  output  at  all  corner  grid  points.  The  model  used  the 
MSC/NASTRAN  "topological"  option  for  grid  point  stress 
calculation  [14].  This  method  assumes  stresses  are  continuous 
across  connecting  elements.  Following  the  stress  intensity 
factor  calculation  technique  described  in  the  computer 
implementation  section,  only  those  stress  grid  points  along  the 
line  of  the  crack,  and  at  a  distance  of  five  to  ten  percent  of 
the  cracks'  length  ahead  of  the  crack  were  used  in  the  stress 
intensity  factor  determination.  The  baseline  model  stress  grid 
points  in  the  crack  tip  area  were  only  .001  inches  apart 
allowing  five  grid  points  in  the  calculations.  The  courser 
finite  element  model  (8,610  DOF)  had  only  two  stress  points  in 
the  calculation  zone,  while  the  fine  model  (15,842  DOF)  had 
nine  grids  in  the  calculation  of  K^. 
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13858  DOF  MASTRflM  MODEL 


Figure  9.  Baseline  MSC/NASTRAN  Finite  Element  Model 


35 


B .  Boundary  Element  Method 

The  boundary  element  model  constructed  consisted  of 
287  elements  resulting  in  574  degrees  of  freedom.  The  model 
is  pictured  in  Figure  10.  The  crack  is  descretized  with  F.R. 
Harris's  refinement  technique  [8]  resulting  in  elements  at 
the  crack  tip  with  a  length  of  0.0005  inches.  The  model  is 
restrained  from  rigid  body  movement  by  fixing  both 
displacement  boundary  conditions  for  the  far  right  element  on 
the  line  of  symmetry.  The  boundary  element  model  had  five 
element  midpoints  in  the  allowable  zone  for  calculation. 

C.  Comparisons 

Stress  Intensity  Factor  calculations  were  completed 
on  all  three  finite  element  models  and  the  boundary  element 
model  using  the  stress  extrapolation  method 

Kj.  =  Limj,^Q[Cy  (2rtr)^/^l  (43) 

2 

The  values  of  were  plotted  against  r  and  r  to  graphically 

determine  Kj  at  r-0.  Linear  regression  fits  were  made  for 

both  fits.  As  was  discussed  in  the  computer  implementation 

2 

section  (Appendix  A) ,  the  r  method  was  necessary  for  cases 
where  the  crack  length  approaches  the  second  hole  as  the  r 
method  yields  poor  curve  fits.  For  this  case,  it  was  not 
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deemed  necessary  as  the  crack  length  is  relatively  short,  but 

was  still  done  for  comparison  purposes.  Final  calculations 

2 

for  Kj  for  the  r  fit  and  the  r  fit  are  given  in  Table  I.  It 
can  be  seen  that  both  fits  gave  essentially  the  same  answers, 
as  was  expected  for  this  case.  Also,  the  boundary  element 
prediction  was  within  two  to  three  percent  of  the  baseline 
finite  element  predictions. 

The  agreement  between  the  boundary  element  model  and 
the  baseline  finite  element  model  is  encouraging  considering 
the  difference  in  degrees  of  freedom  (13858  to  574).  From 
this  simple  statement  the  reader  would  conclude  that  the 
boundary  element  method  is  24  times  more  efficient.  But  the 
user  must  remember  that  the  boundary  element  model  was  a 
"full"  matrix  without  the  banded  symmetry  common  to  the 
finite  element  method.  A  highly  optimized  finite  element 
code,  such  as  MSC/NASTRAN,  has  a  built  in  nodal  resequencer 
to  optimize  the  stiffness  matrix  automatically.  The  VAX 
computer  operates  with  a  "virtual  memory"  scheme.  Matrix 
storage  is  handled  by  writing  to  scratch  files  that  are 
erased  upon  program  completion.  This  makes  it  difficult  to 
compare  storage  requirements  for  both  FEM  and  BEM.  Therefore, 
it  is  instructive  to  examine  the  CPU  times  required  to  run 
all  four  models  as  listed 'in  Table  1.  All  of  the  CPU  times 
are  for  a  Digital  VAX  8350  computer.  It  can  be  seen  that  in 
comparing  the  boundary  element  model  to  the  baseline 
MSC/NASTRAN  model,  it  ran  2.2  times  as  long  even  though  the 
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MSC/NASTRAN  model  used  24  times  as  many  degrees  of  freedom. 
The  boundary  element  model  barely  ran  faster  than  the  8,610 
DOF  course  MSC/NASTRAN  model. 

The  reason  for  the  CPU  time  results  lie  in  the 
relationship  between  matrix  size,  fullness  and  the  CPU  time 
to  invert  and  solve  it.  The  MSC/NASTRAN  Handbook  for  Linear 
Static  Analysis  [13]  outlines  a  relationship  between  problem 
size  and  computer  time.  Basically  the  three  elements  of 
computer  time  are;  overhead  cost,  which  is  dependent  on 
problem  type  but  not  on  problem  size;  initial  matrix  set  up 
costs,  which  involve  computation  of  the  influence  or 
stiffness  matrices;  and  finally  results  costs  which  involve 
solving  the  matrices  for  final  computations.  The  results  cost 
are  the  one  that  increases  rapidly  with  an  increase  in 
problem  size.  Reference  (13]  states  that  for  a  finite  element 
model  with  approximately  100  to  200  grids,  all  three  costs 
are  the  same.  It  is  obvious  that  this  study  has  far  more  than 
200  grids,  so  will  be  dominated  by  the  results  costs. 
Reference  [13]  goes  on  to  give  explicit  formulas  for  CPU 
estimation,  but  the  CPU  formulas  are  proportional  to  the 
number  of  degrees  of  freedom  multiplied  by  the  average  (RMS) 
number  of  active  columns  squared.  The  baseline  MSC/NASTRAN 
output  yielded  a  RMS  value  for  active  columns  after 
resequencing  of  approximately  ninety  colums.  A  full  BEM 
matrix  of  574  by  574  has  a  RMS  column  width  of  332. 
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Therefore : 


CPUbem  *  (574  DOF) (332  columns  RMS ) ^-63268576  (44) 

CPUpEM  “  (13858  DOF) (90  columns  RMS ) ^-112249800 

The  ratio  of  CPU__„  over  CPU__„  is  1.7.  This  indicates  that  a 
preliminary  comparison  of  the  boundary  element  model  to  the 
baseline  MSC/NASTRAN  model  should  have  predicted  a  run  time 
for  the  MSC/NASTRAN  model  of  1.7  times  the  boundary  element 
model,  not  24  times.  (Actual  CPU  time  ratio  was  2.2)  Indeed, 
the  cost  of  a  fully  populated  matrix  is  very  high. 
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Table  I.  MSC/NASTRAN  VS  BEM  Tension  Strip  Results 


Model 

DOF 

Kj  (r  fit) 

Kj  (r^  fit) 

CPU 

(KSI(in)^/^) 

(KSl(in)^/^) 

(min 

Coarse  FEN 

8610 

37.9 

36.3 

96 

Baseline  FEM 

13858 

39.6 

38.9 

183 

Fine  FEM 

15842 

39.5 

38.9 

255 

BEM 

574 

40.4 

40.1 

83 
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IV.  Boundary  Element  vs  p-Verslon  Finite  Element 


PROBE  is  a  commercial  finite  element  code  sold  and 
promoted  by  Noetic  Technologies.  The  code  was  first  conceived 
and  implemented  at  Washington  University's  Center  for 
Computational  Mechanics  in  St.  Louis  under  Dr.  Barna  A. 

Szabo.  The  theoretical  aspects  of  the  p-Version  of  finite 
elements  are  explained  by  Babuska,  Szabo  and  Katz  in 
reference  [21.  The  implementation  of  the  p-Version  into  PROBE 
is  given  by  Szabo  in  reference  [23].  The  innovative  aspect  of 
PROBE  is  that  it  boasts  elements  based  on  variable  order 
polynomials.  By  doing  this,  the  user  can  create  very  rough 
grids  in  the  creation  of  finite  element  models.  By  varying 
the  polynomial  order,  or  p,  increased  accuracy  in  the  results 
is  obtained.  The  second  advantage  of  PROBE  is  that  by  running 
multiple  p  levels  for  a  given  model,  the  user  is  given  an 
indication  of  solution  convergence. 

Noetic  Technologies  worked  with  the  Fort  Worth 
Division  of  General  Dynamics  (GD)  on  a  research  grant  to 
study  the  application  of  the  p-Version  to  a  stress  intensity 
factor  analysis,  and  compare  it  to  a  classical  finite  element 
solution.  A  two  hole  tension  strip,  as  shown  in  Figure  11, 
was  analyzed  by  GD  with  conventional  finite  element  analysis. 
Noetic  Technologies  analyzed  the  same  problem  with  the 
p-Version  PROBE  code  and  published  the  results  in  reference 
[25].  The  GD  model,  shown  in  Figure  12,  involved 
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approximately  1500  degrees  of  freedom.  The  corresponding 

PROBE  model,  as  shown  in  Figure  13,  has  only  29  nodes. 

However,  by  varying  the  value  of  p  from  1  to  8,  the  PROBE 

degrees  of  freedom  varies  between  58  and  1623. 

A  boundary  element  model  was  constructed  of  the 

problem,  as  shown  in  Figure  14,  consisting  of  220  elements  or 

440  degrees  of  freedom.  The  GD  model  predicted  a  of  43.4 

KSI(in)^'^^.  The  PROBE  results  for  p-1  to  p-8  were  plotted  by 

Kj  versus  1/DOF  on  a  semi-logarithmic  scale,  and  the 

resulting  straight  line  extrapolated  to  predict  at  p-«. 

The  final  PROBE  prediction  of  at  p-*  is  43.1  KSI(in)^'^^. 

The  final  BEN  prediction  based  on  a  regression  fit  on  r  was 
1/2 

42.2  KSI(in)  '  .  The  final  BEN  prediction  for  a  regression 

2  1/2 

fit  on  r  was  43.2  KSI(in)  '  .  Both  BBM  predictions  were 

2 

close  to  the  GD  and  PROBE  predictions,  but  the  r  fit  was 
better . 

It  should  be  noted  that  the  PROBE  analysis  gave  an 

indication  of  convergence  to  the  final  answer.  The  BEM  model 
2 

with  an  r  regression  fit  was  almost  exact  in  its' 
correlation  with  the  PROBE  results. 
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width  =  5.0  in 
helght=10  In 
dia  =  1.0  In 
S  =  20  KSI 
pitch  =  2.1in 
a  =  .55  In 


Figure  11.  PROBE  Two  Hole  Tension  Strip  Problem 
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Approximately  1500  degrees  of  freedom 


Figure  12.  GD  Model  of  PROBE  Problem 
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Approximately  29  nodes 


e  13.  PROBE  model  of  PROBE 


V.  Boundary  Element  vs  Bowie  Solution 


Bowie  [3]  studied  the  problem  of  a  crack  growing  from 
a  circular  hole  in  an  infinite  plate  as  shown  in  Figure  15. 
His  solution  is  well  published  and  can  be  shown  in  the  form, 

K  •=  a(  na)^'^^(3  (45) 

where  fi  -  /(a/r).  Other  individuals,  specifically  Grandt, 
Brussat  and  Newman  [1],  have  employed  various  technique  to 
improve  on  Bowies  fi  term.  For  the  example  problem,  <t-46  KSl, 
r>-.125  in,  and  a/r-.5  .  Using  the  value  of  a/r-0.5,  Bowie, 
Brandt,  Brussat  and  Newman  calculate  a  value  of  1.73,  1.735, 
1.733  and  1.728  for  0  respectively  (II .  When  inserted  into 
Eqn  (45),  this  results  in  calculations  of  35.26,  35.26, 
35.32  and  35.22  KSI(in)^/^. 

A  boundary  element  model  was  created  comprising  of 
72  elements.  To  model  an  infinite  domain,  a  different 
modeling  technique  is  required  than  for  the  finite  domains. 
The  model  is  shown  in  Figure  16.  The  model  is  again  a 
representation  of  the  "upper"  half  of  the  geometric  boundary. 
As  described  in  the  Computer  Implementation  section,  a  line 
of  symmetry  is  assumed  along  the  x  axis.  Phantom  "image" 
elements  are  calculated  by  the  TWOFS99  program  for  the  lower 
half.  One  problem  is  the  crack  itself.  Unlike  the  finite 
domain  problems,  the  crack  elements  cannot  be  on  the  line  of 
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symmetry  as  the  program  could  not  distinguish  the  actual 
elements  representing  the  "upper"  face  of  the  crack,  from  the 
"image"  elements  representing  the  "lower"  face  of  the  crack. 
Therefore  the  line  of  elements  representing  the  "upper"  face 
of  the  crack  are  modeled  with  a  small  crack  opening 
offset  as  shown  in  Figure  17.  The  elements  along  the  crack 
are  arranged  in  a  straight  line  between  the  crack  opening 
offset  and  the  crack  tip.  The  "image"  elements  are  therefore 
calculated  with  an  equal,  but  opposite,  location  below  the 
y-0  line  of  symmetry.  The  objective  is  to  model  the  crack 
opening  offset  as  small  as  possible  to  best  represent  the 
actual  crack,  which  has  no  such  offset.  But  the  offset  must 
be  large  enough  for  the  TWOFS99  program  to  differentiate 
between  the  two  faces  of  the  crack.  This  is  usually  a 
function  of  the  accuracy  of  the  computer  the  program  is 
running  on.  The  Bowie  model  uses  an  offset  of  5.0(10)”® 
inches  from  the  y-0  line  of  symmetry,  to  the  intersection  of 
the  "upper"  face  of  the  crack  with  the  circumference  of  the 
hole.  This  results  in  an  initial  offset  five  orders  of 
magnitude  smaller  than  the  actual  y  displacement  at  that 
point.  The  crack  itself  is  again  modeled  with  the  F.R.  Harris 
refinement  technique  which  concentrates  25  elements  in  the 
crack  tip  area.  The  model  is  symmetric  about  the  y-0  axis  by 
imposed  symmetry,  but  the  element  along  the  circumference  of 
the  hole,  opposite  from  the  crack,  is  restrained  from  x 
displacements  to  prevent  rigid  body  translation. 
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Figure  1 


line  of  symmetry 


r  =  radius 

a  =  crack  length 

0  =  crack  opening  (exaggerated) 


.  Infinite  Domain  BEM  Craclt  Modeling  Technique 
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since  there  are  no  elements  in  the  area  o£  the  stress 

field  used  for  calculations,  points  of  data  calculations 

must  be  placed  there.  As  can  be  seen  in  Figure  16,  eleven 

data  points  were  placed  in  the  line  of  the  crack,  at  a 

distance  of  1.05a  to  1.10a.  The  y  stresses  were  recovered  at 

these  points,  and  used  to  create  stress  extrapolation 

predictions  for  K^.  As  before,  both  regression  fits  on  r  and 
2 

r  were  completed. 

Based  on  the  72  element  model,  the  Kj  prediction 

1/2 

based  on  a  regression  fit  on  r  is  35.0  KSI(in)  '  .  The  same 

1/2  2 

model  predicted  35.6  KSI(in)  '  based  on  a  r  fit.  In  this 

2 

instance  the  r  fit  was  more  accurate  than  the  r  fit,  but  the 
significant  observation  is  that  both  methods  provided  a 
prediction  within  one  percent  of  the  analytical  predictions 
of  Bowie,  Grandt,  Brussat  and  Newman  [1]. 
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VI .  Boundary  Element  vs  Shivakumar  Solution 


The  next  problem  attempted  is  an  extension  of  the  Bowie 

problem  of  Section  VI.  A  second  hole  is  added  to  the  Bowie 

problem  to  simulate  the  two  hole  tension  strip  problem  of 

Sections  IV  and  V,  only  the  domain  is  infinite,  not  finite. 

Shivakumar  and  Foreman  solved  this  problem  [19]  with  a  series 

approach  based  on  the  Muskhelishvili  formulation.  The  solution 

is  incorporated  into  the  NASA  crack  growth  computer  program 

NASA/FLAGRO  [15].  By  selecting  a  far  field  stress  of  46  KSI, 

crack  length  of  0.0625  inches,  hole  diameter  of  0.25  inches  and 

a  hole  separation  of  1.0  inch,  the  analytical  prediction  of 

from  the  NASA/FLAGRO  program  is  36.03  KSI(in)^'^^. 

The  analytical  solution  assumes  a  row  of  holes  in  an 

infinite  plate.  To  properly  model  the  geometry  with  boundary 

elements,  three  holes  were  included  in  the  analysis.  This 

included  one  hole  on  either  side  of  the  flawed  hole.  The 

modeling  techniques  were  identical  to  the  Bowie  solution  model 

in  Section  VI.  The  model  is  depicted  in  Figure  18.  The  model 

consisted  of  148  elements,  or  296  degrees  of  freedom.  The 

2 

stress  data,  as  before,  was  fit  to  both  r  and  r  .  The 

1  /2 

prediction  for  r  was  35.3  KSI(in)  '  while  the  prediction  for 

2  1/2  2 
an  r  fit  was  36.0  KSI(in)  '  .  In  this  case  the  r  fit  more 

closely  approximated  the  analytical  solution.  However,  both 

fits  were  within  two  percent  of  the  analytical  solution  with 
2 

the  r  fit  being  only  0.09  percent  different. 
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detail  B>B 


Figure  18.  Boundary  Element  Model  of  Shivakumar  Problem 
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It  is  important  to  observe  again  that  the  analytical 
solution  assumes  an  infinite  row  of  holes.  Obviously  the 
three  holes  nearest  to  the  crack  dominated  the  solution,  but 
additional  refinement  could  be  achieved  by  including  more  of 
the  remaining  holes. 


VII .  Two  Hole  Tension  Strip  Parametric  Study 


The  final  analytical  task  is  a  parametric  study  for  a 
two  hole  tension  strip  analysis.  The  comparisons  to 
conventional  finite  element  analysis  for  this  configuration 
problem  was  established  with  MSC/NASTRAN  in  Section  III,  and 
Noetic  PROBE  in  Section  IV.  Correlation  of  the  boundary 
element  method  and  modeling  techniques  employed  in  this  study 
were  shown  with  the  comparison  to  the  infinite  domain 
problems  of  the  Bowie  solution  in  Section  V,  and  the 
Shivakumar  solution  in  Section  VI.  This  section  is  an 
analysis  of  a  two  hole  tension  strip  with  the  geometry  and 
boundary  conditions  as  shown  in  Figure  19.  The  edge  distance 
from  the  center  of  the  holes  to  the  side,  top  and  bottom 
edges  is  established  as  three  hole  diameters.  All  cases  will 
be  analyzed  for  a  far  field  tension  stress  of  46  KSI.  The 
parameters  that  are  varied  in^this  study  are  hole  diameter, 
d,  hole  separation  (center  to  center)  p  (expressed  as  a  ratio 
of  hole  diameter),  and  crack  length  a.  This  study  expressed 
crack  length  as  a  ratio  where 

crack  ratio  »  a/(  P  -  D  )  (46) 

where  a  -  crack  length  (in) 

D  -  hole  diameter  (in) 

P  -  hole  pitch  as  a  ratio  of  D(  D  in) 
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dia  =  .25,  .33,  .50  In 
pitch  =  3,  4,  5  diameters 
e/d  =  3  dia 

Crack  Ratio  =  a/(pltch  -  dia) 

=  .1,  .2,  .3,  .4,  .5,  .6,  .7,  .8,  .9 
S  =  46  KSI 


Figure  19.  Two  Hole  Tension  Strip  Parametric  Study 


This  enables  the  crack  length  to  be  expressed  as  a  fraction  of 
the  distance  of  material  available  between  the  two  holes. 
Therefore  a  crack  ratio  of  zero  corresponds  to  no  crack  at  all, 
and  a  crack  ratio  of  one  implies  the  crack  has  broken  through 
from  the  first  hole  into  the  second  hole. 

The  study  included  hole  diameters  of  0.25,  0.33  and 
0.50  inches.  The  pitch  was  analyzed  for  3D,  4d  and  5D,  and  the 
crack  ratio  was  analyzed  for  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7, 
0.8,  and  0.9  .  The  values  for  hole  diameter,  edge  distance  and 
hole  separation  were  chosen  to  represent  realistic  geometry 
found  in  actual  applications.  The  final  study  involved  81 
models  of  the  different  configurations  listed  here,  as  well  as 
18  more  models  for  additional  work  not  included  in  the  baseline 
analysis . 

The  size  of  the  models  varied  from  240  to  340  elements. 
All  of  the  models  were  created  by  the  same  model  generator, 
CHOLE,  as  documented  in  Appendix  A.  The  crack  tip  refinement 
method  was  the  F.R.  Harris  technique  (8J. 

The  stress  field  data  was  collected  at  a  distance  five 
to  ten  percent  of  the  crack  length  ahead  of  the  crack  tip.  This 
is  the  same  method  used  throughout  this  thesis.  The  stress 
field  is  used  to  predict  the  mode  I  crack  tip  stress  intensity 
factor,  Kj,  by  using  the  equation 

Kj  -  Limj._^Q  (ay  (2nr)^/^J 


(43) 


as  documented  in  Appendix  A  for  the  program  TWOFS99_EX. 

TWOFS99_EX  extracted  the  stress  field  data  for  all  of  the 

boundary  element  models,  computed  the  values  of  and  r,  and 

2 

fit  the  data  with  a  linear  regression  analysis  of  vs  r  . 

Throughout  this  thesis,  predictions  based  on  regression  fits 
2 

of  r  and  r  have  been  presented.  The  results  were  for  the 

geometry  analyzed.  There  was  no  significant  difference  in  which 

fit  was  chosen,  and  neither  regression  fit  was  consistently 

more  accurate  than  the  other.  However,  during  the  course  of 

this  study,  it  was  found  that  for  crack  ratios  approaching  0.9, 

due  to  the  influence  of  the  approaching  second  hole,  the  vs 

r  curve  is  decidedly  non-linear.  Therefore,  a  linear  regression 

2 

fit  was  non-representative.  The  vs  r  curve  was  much  more 

linear,  and  the  regression  fit  of  that  data  was  representative. 

For  this  reason,  all  predictions  presented  in  this  section 

2 

are  based  on  a  vs  r  regression  fit  only.  Examples  of 

2 

data  plotted  against  r  and  r  are  presented  in  Appendix  F. 

The  Kj  calculations  for  hole  diameters  of  0.25,  0.33 
and  0.50  inches  are  presented  in  Figures  20,  21,  and  22 
respectively.  The  calculated  Kj  values  are  plotted  against  the 
crack  ratios  and  are  presented  as  a  family  of  curves  varying  by 
the  pitch.  The  figures  show  the  trend  is  for  increasing  values 
of  Kj  for  increasing  crack  ratios,  and  for  increasing  for 
increasing  hole  diameter.  The  also  increased  for  increasing 
pitch  ratios.  If  the  data  were  plotted  on  the  same  graph,  the 
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(KSi  SQR(IN) 


STRESS  INTENSITY  FRCTQR  VS  CRHCK  RATIO 


FOR  HOLE  DIAMETER  =  0.33  INCHES 


CRACK  RATIO  (A/CPITCH  -  DIA]  ) 


Figure  21.  vs  Crack  Ratio  for  Hole  Diameter  -  0.33 
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(KSl  SQR(IN) 


STRESS  INTENSITY  FRCTOR  VS  CRACK  RATIO 


FOR  HOLE  DIRriETER  =  0-50  INCHES 


Figure  22.  vs  Crack  Ratio  for  Hole  Diameter  ■  0.50 
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presentation  would  be  confusing  as  the  curves  would  overlap, 
making  interpretation  of  results  difficult.  A  better  way  of 
presenting  the  data  from  this  study  is  the  Stress  Intensity 
Factor  Correction  Coefficient,  0,  as  defined  by 

0  -  Kj/<T(na)^/^  (46) 

By  "normalizing"  the  stress  intensity  factor,  the  influence 
of  far  field  stress  and  crack  length  are  removed,  allowing 
for  isolation  of  the  geometric  Correction  Coefficient  (0)  of 
the  problem  being  solved.  The  0  factors  are  presented  in 
Figure  23.  It  was  found  that  by  plotting  0  versus  the  crack 
ratio,  a  family  of  curves  varying  by  the  pitch  ratio  could  be 
produced.  Once  plotted  with  these  parameters,  the  variation 
of  0  with  the  hole  diameter  was  found  to  be  invariant.  Thus, 
Figure  23  represents  a  useful  tool  in  the  analysis  of  the  two 
hole  tension  strip  with  the  edge  constraints  presented  in  the 
beginning  of  this  section.  The  values  of  all  computed  0 
factors  for  all  of  the  models  run  are  presented  in  Tables  II, 
II,  and  IV. 

It  is  interesting  to  note  the  compression  of  the  0 
curves  at  the  higher  pitch  ratios,  at  crack  ratios  above  0.5. 
To  analyze  this  phenomena,  for  a  hole  diameter  of  0.25 
inches,  two  tiditional  curves  with  a  pitch  ratio  of  3.5  and 
4.5  were  created.  These  curves  were  plotted  with  the  previous 
0  curves  to  create  Figure  24.  This  shows  that  there  is  a 
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BETR  FHCT0R  VS  CRRCK  RRTI0 


Figure  23.  0  Factor  vs  Crack  Ratio 


KI/  iSIGry  SQRTIPI  h] 


BETn  FRCT0R  VS  CRRCK  RRTI0 

HOLE  DIRMETER  =  0-25  INCHES 


Figure  24.  Factor  vs  Crack  Ratio  for  Hole  Diameter-O . 25 
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compression  in  the  p  curves  at  the  location  mentioned 

earlier.  It  was  postulated  that  this  was  a  "net  area"  effect 

relating  to  the  rigid  edge  distance  criteria  of  the  original 

problem.  For  a  given  diameter  hole,  the  pitch  was  varied  as  a 

ratio  of  the  diameter,  but  the  edge  distances  remained 

constant  at  three  diameters.  Therefore,  as  the  crack  ratio 

grows  towards  0.9,  the  reduction  in  net  area  as  a  percentage 

of  the  total  original  pre-cracked  net  area,  is  higher  for  the 

higher  pitch  ratios.  The  effects  of  this  would  be  increased 

as  the  crack  grew  in  length.  To  examine  this  trend,  the  3 

factors  from  Figure  23  were  modified  to  calculate  (3  based  on 

net  stress,  instead  of  far  field  stress,  a,  and  then 

calculate  as  follows 

net 

»net  -  ■'l/'W 

The  results  are  shown  in  Figure  25.  Both  (5  and  Pj^g^  are 
plotted  against  crack  ratio.  The  plot  shows  that  at  crack 
ratios  above  0.5,  as  the  fi  factors  based  on  far  field  stress 
began  to  increase  uniformly  in  value,  the  Pj^g^  factors  based 
on  net  stress  cross  over  as  the  effects  of  pitch  ratio  seem 
to  reverse.  It  is  further  postulated  by  the  author,  that  if 
the  net  section  effects  were  subtracted  from  the  final  3 
curves  of  Figure  23,  a  family  of  0  curves  would  thus  be 
created  with  the  same  generic  trends  of  Figure  23,  but 
without  the  collapse  of  curves  at  the  higher  pitch  ratios 
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BETA  -  KI/  (SIGY 


BETA  RND  BETA  -NET  FfiCTORS  VS  CRACK  RATIO 


FOR  HOLE  DIAMETER  =  0-25  INCHES 


0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0 

CRACK  RATIO  (A/(PITCH  -  DIA3) 


Figure  25.  (3  and  .  Factors  vs  Crack  Ratio 
^  net 
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above  a  crack  ratio  of  0.5  . 


It  can  be  seen  from  Figure  23  that  for  crack  ratios 
up  to  0.5,  the  effects  of  the  initial  hole  are  dominate,  with 
the  influence  of  the  hole  decreasing  with  increased  distance 
from  the  hole.  At  a  crack  ratio  of  0.5,  the  crack  begins  to 
approach  the  second  hole  and  the  value  of  fi  now  increases 
with  the  decrease  in  distance  to  the  second  hole.  So  all  of 
the  items  of  the  initial  problem  can  be  seen  in  the  final  fi 
curves  of  Figure  23.  The  first  hole  is  seen  in  the  high 
initial  values  of  with  the  effects  of  the  hole  decreasing 
with  distance.  The  (3  value  are  at  a  minimum  approximately 
half  way  between  the  holes,  with  the  effects  of  the  second 
hole  seen  as  the  0  values  Increasing  with  the  crack  tip 
approaching  the  second  hole.  The  effects  of  the  edge 
distances  are  seen  in  the  "collapse"  of  the  0  curves  at  high 
pitch  and  crack  ratios.  All  of  this  is  in  addition  to  the 
obvious  effects  of  pitch  and  crack  ratio  as  a  function  of 
hole  diameter. 


69 


Table  II.  Parametric  Study  p  Factors  for  Pitch  Ratio  -  3  Dia 


Crack  Ratio 

^dia-0.50 

^dia-0.33 

^dia-0 . 

.1 

2.07 

2.07 

2.07 

.2 

1.62 

1.62 

1.62 

.3 

1.43 

1.43 

1.43 

.4 

1.34 

1.34 

1.34 

.5 

1.31 

1.31 

1.31 

.6 

1.31 

1.31 

1.31 

.7 

1.34 

1.35 

1.34 

.8 

1.46 

1.46 

1.48 

.9 

1.90 

1.90 

1.90 

25 
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Table  IV.  Parametric  Study  fi  Factors  for  Pitch  Ratio  -  5  Dia 


Crack  Ratio 

^dia-0.50 

^dia-0.33 

^dia»0 . ; 

.1 

1.55 

1.55 

1.55 

.2 

1.26 

1.25 

1.25 

.3 

1.17 

1.17 

1.17 

.4 

1.14 

1.14 

1.14 

.5 

1.15 

1.15 

1.15 

.6 

1.18 

1.18 

1.18 

.7 

1.23 

1.23 

1.26 

.8 

1.32 

1.32 

1.34 

.9 

1.58 

1.58 

1.58 
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VIII.  Parametric  Study  Application 


This  section's  purpose  is  to  present  an  example  of  how 
the  fracture  mechanics  engineer  in  the  aircraft  industry  might 
apply  the  results  of  the  parametric  study  undertaken  in  the 
last  section.  The  example  outlined  here  is  hypothetical  and  not 
intended  to  limit  the  potential  usage  of  the  stress  intensity 
data  in  the  last  section. 

It  will  be  assumed  that  there  is  a  requirement  for  a 
particular  structure,  in  this  case  a  machined  fitting  made  out 
of  7075-T6  aluminum  plate  (with  E=10300  KSI,  and  v«0.33).  The 
fitting  is  to  have  a  service  life  of  500  flight  hours.  It  will 
be  further  assumed  that  the  only  significant  load  on  the 
fitting  is  aircraft  pressurization,  and  therefore  the  fitting 
will  experience  one  load  cycle  per  flight.  The  average  flight 
for  this  airplane  will  be  one  hour. 

To  establish  the  Damage  Tolerance  of  this  part  under 
current  Air  Force  requirements  [7],  this  analysis  will  qualify 
the  fitting  as  being  "slow  crack  growth"  structure,  and 
therefore  it  must  be  shown  that  two  service  lifetimes  of  slow 
crack  growth  exist. 

Though  fracture  mechanics  and  stress  intensity  factor 
calculations  are  based  on  theory,  fatigue  crack  growth  is 
empirical.  Crack  growth  for  a  particular  stress  cycle  is  a 
function  of  the  change  in  stress  intensity  and  stress  ratio 
(the  minimum  stress  divided  by  the  maximum  stress  in  a  cycle) 
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for  a  given  material.  This  thesis  will  not  cover  the  theory  of 
crack  growth  analysis,  nor  of  material  fracture  properties. 
However,  it  is  important  for  the  reader  to  understand  that  for 
a  given  material,  the  crack  growth  increment  for  a  given  stress 
cycle  is  dependent  on  the  stress  intensity  at  the  time  of  load 
application. 

Many  software  codes  have  been  written  to  do  fatigue 
crack  growth  analysis  (CRACKS,  CRKGRO,  FLAGRO,  etc.),  and  they 
all  share  certain  traits  in  common.  After  input  of  basic 
material  fracture  properties  for  the  material  being  used,  the 
stress  spectrum  is  input.  Then  the  algorithm  to  calculate  the 
stress  intensity  factor  throughout  the  analysis  is  selected. 
(The  stress  intensity  factor  will  vary  with  the  crack  length 
and  applied  stress)  Most  crack  growth  codes  have  a  library  of 
predefined  crack  stress  intensity  solutions  to  choose  from. 

Most  codes  also  allow  the  user  to  input  a  "look-up"  table  of 
stress  intensity  data  vs  crack  length.  The  look-up  data  is 
usually  in  the  form  of  a  (3  factor,  as  calculated  in  the  last 
section.  The  format  of  the  stress  spectrum  is  usually  written 
as 


"max'  "min' 


(48) 
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where 


»  maximum  stress 

max 

~  “  minimum  stress 

min 

cycles  *  number  of  repetitions 

The  information  of  equation  (48)  can  be  repeated  to  create 
layers  in  a  complex  stress  spectrum.  The  spectrum  applied  in 
this  analysis  is  very  simple  as  it  has  only  one  cycle  per 
flight.  Stress  analysis  of  the  fitting  indicated  an  applied 
stress  of  30  KSI  under  fully  pressurized  conditions,  with  0  KSI 
unpressurized.  Therefore  the  stress  spectrum  for  one  flight 
would  be 

-  30  KSI  (49) 

-  0  KSI 

-  1 

Most  engineers  attempt  to  compile  a  spectrum  into  a  "bloclc" 

that  would  represent  many  flights,  and  then  repeat  the  block 

until  the  service  life  requirements  are  met.  This  analysis 

defines  100  flights  to  be  a  block,  therefore  one  block 

represents  100  flight  hours  of  life.  Two  service  lives  of  slow 

crack  growth  must  be  shown  before  critical  crack  length  is 

reached.  Critical  crack  length  is  either  loss  of  a  part,  or 

when  the  crack  length  grows  to  a  point  where  the  local  stress 

intensity  factor  for  a _ exceeds  the  material  fracture 

max 


or 

max 

‘’^rain 

cycles 
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toughness.  (This  thesis  will  also  not  cover  the  Air  Force 
residual  strength  requirements)  To  achieve  two  lifetimes  of 
slow  craclt  growth,  an  assumed  initial  crack  must  not  grow  to 
critical  crack  length  before  1000  flight  hours. 

In  this  example,  the  fatigue  crack  growth  computer 
program  NASA/FLAGRO  [15]  was  used.  This  was  also  the  source  for 
the  calculations  of  the  Shivakumar  solution  used  in  section 
VII.  The  built  in  material  fracture  data  for  7075-T6  aluminum, 
and  a  constant  spectrum  of  0  -  30  KSI  was  used  for  all  versions 
of  this  analysis. 

Three  approaches  were  taken  in  the  analysis  of  the 
fitting.  It  was  assumed  the  critical  crack  location  was  a 
through  the  thickness  flaw  emerging  from  a  fastener  hole,  with 
a  geometry  as  shown  in  Figure  19.  The  fitting  has  0.25  inch 
fastener  holes  with  a  hole  separation  (pitch)  of  four  diameters 
(1.0  inch  in  this  case).  The  initial  flaw  sizes  are  dictated  by 
the  Air  Force,  and  vary  by  type  and  location.  The  size  is 
determined  by  the  largest  "rogue  flaw"  that  could  be  induced  in 
the  fitting  during  manufacture,  assembly,  or  service  use  that 
could  not  be  detected  by  routine  non-destructive  inspections 
(NDI)  with  a  90  percent  probability  of  detection,  and  a  95 
percent  confidence.  It  was  assumed  here  that  the  local  NDI  was 
not  very  good,  and  that  an  initial  through  the  thickness  crack 
size  of  0.075  inches  would  be  used.  This  is  convenient  as  this 
translates  into  a  crack  ratio  of  0.1  (using  the  definitions  of 
the  previous  section). 
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The  fitting  was  analyzed  using  three  different 
approaches  to  the  calculation  of  the  stress  intensity  factor  as 
a  function  of  crack  length  and  applied  stresses.  The  first 
method  used  the  Bowie  solution  approach  to  idealize  the  fitting 
as  a  hole  in  an  infinite  plate.  The  spectrum  was  applied  to  the 
initial  flaw  and  grown  to  a  length  of  0.75  inches  which 
represents  the  length  required  to  "break  through"  into  the 
second  hole.  The  second  method  used  the  Shivakumar  solution 
assuming  a  row  of  holes  in  an  infinite  plate.  This  analysis 
also  terminated  upon  the  crack  reaching  the  second  hole.  The 
last  method  involved  the  3  factors  derived  in  the  last  section. 
The  3  factors  were  placed  in  a  3  look-up  table  as  a  function  of 
crack  length.  Stress  intensity  factors  were  then  calculated  for 
a  given  crack  length,  a,  and  a  given  applied  stress,  a,  as 
follows 

Kj  -  0  (na)^/^  3  (50) 

The  results  of  the  three  analysis  are  shown  in  Figure  26.  It 
can  be  seen  that  the  Bowie  solution  method  was  the  least 
conservative,  as  it  did  not  consider  the  second  hole,  or  the 
tension  strip  edge  effects.  The  Shivakumar  solution  method  was 
the  second  least  conservative  as  it  did  not  consider  the  finite 
edge  effects.  Both  the  Bowie  and  Shivakumar  methods  grew  the 
crack  until  it  reached  the  second  hole.  The  last  method,  or  "3 
look-up"  table  method  was  the  most  conservative.  The  crack  did 
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CRACK  LENGTH  A  ( fN) 


FRTIGUE  CROCK  CURVES  F0R  FITTING  RNRLYSIS 


Figure  26.  Example  Analysis  Fatigue  Crack  Growth  Curves 
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not  reach  the  second  hole  as  the  local  crack  tip  stress 
intensity  factor  reached  the  material  fracture  toughness  at  a 
crack  length  of  only  0.66  inches.  The  look-up  method  also  was 
the  only  one  unable  to  show  1000  hours  of  slow  crack  growth, 
thus  not  meeting  the  design  requirements. 

From  this  simple  example,  it  can  be  seen  how  detailed 
analysis  through  the  fi  look-up  table  method  enables  an  engineer 
to  analyze  detailed  geometry  beyond  the  scope  of  the  common 
solutions  found  in  most  fatigue  crack  codes.  In  this  case,  it 
would  have  been  unconservative  to  ignore  the  effects  of  the 
second  hole,  and  the  edge  effects.  The  output  from  the 
NASA/FLAGRO  program  are  included  in  Appendix  E. 
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IX.  Conclusions 


Application  of  the  Boundary  Element  Method  to 
structures  problems  is  just  beginning  to  become  popular  in 
the  aircraft  industry.  Traditional  Finite  Element  Methods  are 
still  the  predominate  technique  used.  However,  the  Finite 
Element  Method  is  expensive  in  both  manpower  and  computer 
costs,  and  cost  saving  alternatives  are  always  being  sought. 

The  Fictitious  Stress  Method,  presented  in  this 
thesis,  is  shown  to  correlate  well  with  both  analytical  and 
FEM  solutions.  The  BEM  was  shown  to  work  well  for  the 
parametric  study  of  Section  VII.  A  complicated  fracture 
mechanics  problem  with  no  analytical  solution  was  solved  for 
various  geometry,  with  the  results  displayed  as  a  family  of 
0-curves  in  Figure  23.  These  curves  in  themselves  are 
important  as  they  represent  useful  Stress  Intensity  Factor 
correction  factors  for  the  various  geometric  configurations 
analyzed  in  Section  VII. 

It  has  been  shown  herein  that  BEM  is  an  acceptable 
method  for  fracture  mechanics  analysis  and  can  be  used  in 
fatigue  crack  growth  predictions  for  Air  Force  Durability  and 
Damage  Tolerance  Analysis  (DADTA).  The  analysis  in  Section 
VIII  shows  how  easily  the  results  of  the  BEM  work  in  Section 
VII  could  be  applied  to  a  "real"  design  problem  and  prevent 
unconservative  structural  life  predictions. 
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One  possible  source  for  additional  work  is  in  the 
area  of  FEM  and  BEN  combined  in  a  single  solution.  This  might 
prove  to  be  the  best  of  both  worlds  with  a  fine  grid  FEM 
model  near  the  crack  tip,  and  a  coarse  BEN  definition  of  the 
external  boundaries  of  a  problem.  Additional  study  should  be 
done  to  see  if  the  BEM/FEM  combined  analysis  offers 
advantages  in  actual  applications  to  each  method  used 
separately. 

The  question  of  increased  efficiency  is  a  more 
difficult  one.  Though  dramatic  reductions  in  degrees  of 
freedom  are  shown  for  comparable  accuracy  of  analysis,  final 
CPU  time  is  not  always  improved.  Since  the  CPU  time  is 
proportional  to  the  square  of  the  Root  Mean  Square  (RMS) 
number  of  active  columns  multiplied  by  the  total  degrees  of 
freedom  (DOF),  the  BEM  would  have  a  CPU  time  advantage  for 
smaller  problems  where  the  DOF  factor  would  dominate  the 
squared  RMS  term.  This  indicates  that  the  BEM  is 
computationally  more  efficient  for  problems  up  to  a  certain 
size.  Even  at  problem  sizes  of  13858  DOF  the  BEN  still  has  a 
CPU  time  advantage  of  1.7  (reference  Section  III). 

It  should  also  be  noted  that  the  BEN  program  used  in 
this  thesis  utilized  only  single  precision  accuracy  which  on 
the  VAX  computer  provides  six  significant  digits.  This  helped 
improve  the  BEM  computer  efficiency  and  still  obtain  the 
excellent  correlation  to  the  FEM  and  analytical  results 
documented  in  this  thesis. 


It  is  therefore  concluded  that  for  the  structural 
fracture  mechanics  problems  analyzed  in  this  thesis,  the  BEM 
accurately  derived  Stress  Intensity  Factors  for  fracture 
analysis,  and  produced  a  minimum  computational  efficiency 
improvement  of  1.7  over  traditional  FEM. 
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Appendix  A:  Computer  Implementation 


The  source  for  the  computer  program  used  in  this 
study  was  a  FORTRAN  program,  TWOFS,  for  the  fictitious  stress 
boundary  element  method  published  by  Crouch  and  Starfield 
[4].  The  version  used  in  this  study  was  converted  to  the 
Microsoft  BASIC  computer  language  for  ease  of  implementation 
on  PC  class  computers.  Upon  initiation  of  actual 
calculations,  it  was  decided  to  port  the  program  up  to  a 
Digital  VAX  computer  for  speed  purposes,  so  limited  code 
changes  were  made  to  run  under  VAX  BASIC  3.1.  Additional 
small  changes  were  made  to  facilitate  post  processing  by 
outputting  desired  calculations  to  an  external  file.  The  VAX 
Basic  version  of  TWOFS  was  labeled  TWOFS99.  The  final  work 
was  done  on  a  Digital  VAX  8800  running  the  VMS  (V4.7) 
operating  system.  The  average  BEM  model  consisting  of  300 
elements  (600  by  600  matrix  of  influence  coefficients)  took 
approximately  25  minutes  of  CPU  time.  The  CPU  comparisons 
made  to  the  MSC/NASTRAN  finite  element  program  were  done  on  a 
Digital  VAX  8350  computer  as  it  was  the  only  machine  set  up 
to  run  NASTRAN. 

The  flow  of  the  program  TWOFS99  is  identical  to  the 
original  FORTRAN  TWOFS  code.  The  sizes  of  the  matrices  were 
increased  to  allow  larger  problems.  The  TWOFS99  input  was 
modified  to  allow  for  the  BEM  model  to  be  read  from  a  disk 
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I 

file.  This  was  particularly  important  as  the  final  tension 
^  strip  parametric  analysis  required  over  one  hundred  models  to 

be  built,  run  and  analyzed.  The  models  were  constructed  by  an 
independent  VAX  Basic  program,  CHOLE,  and  written  to  disk  in 
the  format  required  by  TWOFS99.  Aside  from  the  normal  TWOFS 
output,  which  was  also  written  to  a  disk  file,  a  third  file 
was  created  by  TWOFS99  of  unlabeled  final  stress  results.  A 
third  VAX  Basic  computer  program,  TWOFS99_EX,  extracted 
necessary  stress  data  from  the  post  processing  file  created 
by  TWOFS99  and  computed  a  value  for  the  stress  intensity 
factor  based  on  a  regression  fit  technique.  This  allowed  for 
a  great  degree  of  mechanization  in  the  analysis  process. 

A.  Fictitious  Stress  Method  Program  (TWOFS99) 

The  input  file  for  TWOFS99  defines  the  geometry  of 
the  problem,  along  with  the  necessary  boundary  conditions. 

The  program  first  reads  in  values  for  NUMBS,  NUNOS,  KSYM,  PR 
and  E.  NUMBS  defines  the  number  of  straight  line  segments 
which  will  be  input.  NUMOS  defines  the  number  of  additional 
segments  to  establish  data  points  for  displacement  and  stress 
calculations  within  the  body  to  be  analyzed.  KSYM  is  a  code 
to  take  advantage  of  any  lines  of  symmetry  in  a  model  by 
calculating  image  elements  as  mirrored  across  the  line  of 
symmetry  so  that  their  effects  are  included  in  the  final 
results.  KSYM  equal  to  one  implies  no  symmetry  exists,  which 
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was  used  primarily  in  this  study.  KSYM  equal  to  two  implies 
symmetry  about  the  y  axis  at  a  line  x»XSYM.  KSYM  equal  to 
three  implies  symmetry  about  the  x  axis  at  a  line  y»YSYM. 

And  KSYM  equal  to  four  implies  two  axis  of  symmetry  about 
x=XSYM,  y«YSYM.  If  symmetry  is  requested,  the  value  of  XSYM 
and,  or  YSYM  is  input.  PR  is  the  Poison's  Ratio  and  E  is  the 
Young's  Modulus  for  the  material  for  the  problem.  The  field 
stresses  are  next  input  as  PXX,  PYY  and  PXY.  All  input  must 
be  in  consistent  units.  All  input  is  echoed  in  the  output 
f  ile . 

At  line  460  in  the  code,  a  loop  is  entered  from  1  to 
NUMBS.  For  each  iteration  of  the  loop,  values  for  2NUM,  XBEG, 
YBEG,  XEND,  YEND,  KODE,  BVS  and  BVN  are  input.  XBEG,  YBEG, 
XEND  and  YEND  define  the  x  and  y  co-ordinates  for  the 
beginning  and  end  of  the  current  line  segment.  ZNUM 
subdivides  the  current  line  segment  into  that  many  equal 
length  boundary  element  segments.  BVS  and  BVN  are  the 
boundary  conditions  for  all  of  the  boundary  elements  defined 
for  the  current  line  segment,  in  the  shear  and  normal  local 
co-ordinates  of  the  elements  respectively.  KODE  defines  if 
BVS  and  BVN  are  displacement  or  stress  boundary  conditions. 
Remember,  it  is  allowable  to  mix  them  as  indicated  in  [4]. 
KODE  equal  to  one  means  both  are  stresses,  two  means  both  are 
displacements,  three  means  a  shear  displacement  with  a  normal 
stress,  and  four  is  a  shear  stress  with  a  normal 
displacement.  Upon  completion  of  the  loop,  all  input  of  the 
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data  for  the  definition  of  the  boundary  elements  and  boundary 
conditions  is  completed. 

At  line  690  in  the  code,  a  similar  loop  is  entered 
from  1  to  NUMOS .  Here  the  variables  EXTERNL(N, i ) ,  i*l  to  4, 
and  NUMMTX(N)  are  read.  The  data  for  the  interior  points  are 
stored  in  matrices  to  facilitate  changes  made  to  output 
formats.  In  order, the  XBEG,  YBEG,  XEND,  YEND,  and  NUMPD  are 
input  and  placed  in  the  EXTERNL  and  NUMMTX  arrays.  XBEG, 

YBEG,  XEND  and  YEND  are  as  for  the  boundary  element  line 
segment  definitions.  NUMPD  defines  the  number  of  straight 
equally  spaced  points  between  and  including  XBEG,  YBEG,  XEND, 
YEND  to  be  included  for  displacement  and  stress  calculations 
after  the  fictitious  stresses  are  solved  for. 

Lines  1250  through  2000  make  various  calls  to 
subroutines  to  calculate  the  influence  coefficients  for  all 
of  the  boundary  elements,  and  assembles  them  into  a  matrix  C. 
Line  2020  calls  a  Gauss  Elimination  subroutine  to  solve  for 
the  fictitious  stresses  which  are  stored  in  the  matrix  P. 

Line  2100  enters  a  loop  to  calculate  the  unknown  boundary 
conditions  at  all  of  the  boundary  element  midpoints.  And, 
finally,  line  3060  is  a  loop  to  calculate  influence 
coefficients  and  the  resulting  stresses  and  displacements  at 
all  of  the  interior  data  points. 

Line  2920  begins  a  loop  to  store  all  stresses 
computed  at  boundary  element  midpoints,  along  with  the  x 
value  of  the  element  midpoint.  This  data  is  written  to  a  disk 
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file  for  post  processing  in  line  2965.  The  program  TWOFS99  is 
listed  in  Appendix  B. 


B .  Boundary  Element  Generation  (CHOLE) 


This  program  was  written  specifically  for  the  two 
hole  tension  strip  analysis.  The  boundary  conditions  for  the 
study  were  incorporated  into  the  program.  The  user  inputs  a 
problem  title,  hole  diameter,  hole  spacing  and  crack  length. 
The  program  divides  the  crack  length  into  boundary  elements 
using  the  F.R.  Harris  refinement  technique  [8],  and  then 
creates  elements  to  model  the  remainder  of  the  tension  strip 
boundary.  The  final  result  is  a  disk  file  in  the  format 
required  by  TWOFS99  for  analysis.  The  program  CHOLE  is  listed 
in  Appendix  D. 


C .  Stress  Intensity  Factor  Calculation  (TWOFS99  EX) 


The  assumptions  used  are  to  calculate  the 
intensity  factor  for  a  given  problem  by  using  the 
stresses,  a,  normal  to  the  line  of  the  crack.  The 
the  stress  intensity  factor  is  calculated  with  the 
stresses  with  the  equation 


stress 
tension 
value  of 
tension 


K 


I 


a 


( 2nr ) 


1/2 


(42) 
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file  for  post  processing  in  line  2965.  The  program  TWOFS99  is 
listed  in  Appendix  B. 

B.  Boundary  Element  Generation  (CHOLE) 

This  program  was  written  specifically  for  the  two 
hole  tension  strip  analysis.  The  boundary  conditions  for  the 
study  were  incorporated  into  the  program.  The  user  inputs  a 
problem  title,  hole  diameter,  hole  spacing  and  crack  length. 
The  program  divides  the  crack  length  into  boundary  elements 
using  the  F.R.  Harris  refinement  technique  [8],  and  then 
creates  elements  to  model  the  remainder  of  the  tension  strip 
boundary.  The  final  result  is  a  disk  file  in  the  format 
required  by  TWOFS99  for  analysis.  The  program  CHOLE  is  listed 
in  Appendix  D. 

C .  Stress  Intensity  Factor  Calculation  (TWOFS99  EX) 

The  assumptions  used  are  to  calculate  the  stress 
intensity  factor  for  a  given  problem  by  using  the  tension 
stresses,  o,  normal  to  the  line  of  the  crack.  The  value  of 
the  stress  intensity  factor  is  calculated  with  the  tension 
stresses  with  the  equation 

Kj  -  (T  (2nr)^/^  (42) 
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All  of  the  tension  strip  parametric  study  results  were 
processed  through  TWOFS99_EX  for  calculations,  and  the 
results  presented  in  that  section  of  this  report.  The  program 
TWOFS99  EX  is  listed  in  Appendix  C. 
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Appendix  B:  Computer  Program  TWOFS99 


This  appendix  contains  the  listings  of  the  computer 
program  used  for  the  boundary  element  analysis  in  this 
thesis.  The  boundary  element  analysis  program  TWOFS99  was 
basically  extracted  from  Crouch  and  Starfield  (Reference  [4]) 
with  changes  to  output  a  post  processing  file  for  TWOFS99_EX 
to  do  regression  analysis  for  predictions.  The  source 
document  program  was  written  in  FORTRAN,  and  that  was 
converted  into  BASIC.  Also,  the  data  matrix  limits  were 
raised  to  analyze  larger  problems.  The  program  TWOFS99  was 
compiled  under  VAX  BASIC  3.1  . 


10  REM  BOUNDARY  ELEMENT  PROGRAM  TWOFS  5  OCT  85 
20  REM  MODIFIED  FOR  LARGE  MODELS  FOR  VAX 
30  REM 

32  DIM  C(600,600) ,B(600) ,P(600) 

40  DIM  XM(300),YM(300),A(300), 

COSBET(300) ,SINBET(300) ,KOD(300) 

50  DIM  EXTRNL( 300,4) ,NUMMTX( 300) ,OUTPT( 300,10) 

60  REM 

70  REM  PRINT"****************************************" 

80  REM  PRINT"  " 

90  REM  PRINT"  BOUNDARY  ELEMENT  PROGRAM  " 

100  REM  PRINT"  " 

110  REM  print"*****************************************" 
120  REM  PRINT"  " 

130  REM  PRINT"  " 

132  REM  INPUT  "ENTER  INPUT  FILE  NAME  ",QIN$ 

134  REM  INPUT  "ENTER  OUTPUT  FILE  NAME  ",QOUT$ 

136  OPEN  "QIN"  FOR  INPUT  AS  #1 

138  OPEN  "QOUT"  FOR  OUTPUT  AS  #2 

139  OPEN  "QMAT"  FOR  OUTPUT  AS  #3 

140  INPUT  #1,TITLE$ 

150  INPUT  #1, NUMBS, NUMOS,KSYM, PR, E 


160 

IF 

KSYM-1 

THEN 

GOTO 

200 

170 

IF 

KSYM-2 

THEN 

GOTO 

210 

180 

IF 

KSYM-3 

THEN 

GOTO 

240 

190 

IF 

KSYM-4 

THEN 

GOTO 

260 

90 


200 

210 

220 

240 

250 

260 

280 

290 

300 

310 

320 

330 

340 

360 

370 

380 

390 

400 

410 

415 

420 

430 

435 

440 

450 

460 

470 

480 

490 

500 

510 

520 

530 

540 

550 

560 

570 

580 

590 

600 

610 

620 

630 

640 

650 

655 

660 

670 

680 

690 

700 


GOTO  300 
INPUT  #1,XSYM 
GOTO  300 

INPUT  #1,YSYM 
GOTO  300 
INPUT  #1,XSYM 
INPUT  #1,YSYM 
REM 
REN 

INPUT  #1,PXX 
INPUT  #1,PYY 
INPUT  #1,PXY 
REM 

CNST-1.0/(4.0*PI’*(1.0-PR)  ) 

COND-(1.0+PR)/E 

PR1-1.0-2.0*PR 

PR2-2.0*(1.0-PR) 

PR3-3.0-4 .0*PR 
REM 

REM  PRINT"  " 

REM  PRINT"  DEFINE  LOCATIONS,  SIZES,  ORIENTATIONS  AND 
BOUNDARY  CONDITIONS  " 

REM  PRINT"  OF  BOUNDARY  ELEMENTS  " 

REM  PRINT"  " 

REM 

NUMBE-0 

FOR  N-l  TO  NUMBS 

I NPUT  # 1 , 2NUM , XBEG , YBEG , XEND , YEND , KODE , BVS , BVN 
XD- ( XEND-XBEG ) /ZNUM 
YD- ( YEND- YBEG ) /ZNUM 
SW-SQR(XD*XD+YD*YD) 

REM 

FOR  NE-1  TO  ZNUM 
NUMBE-NUMBE+1 
N-NUMBE 

XM( M ) -XBEG+ . 5* ( 2 . 0*NE-1 , 0 ) *XD 
YM ( M ) -YBEG+ . 5  * ( 2 . 0  *NE-1 . 0 ) * YD 
A(M)-. 5*SW 
SINBET(M)-YD/SW 
COSBET(M)-XD/SW 
KOD{M)-KODE 
MN-2*M 
MS-MN-1 
B(MS)-BVS 
B(MN)-BVN 
NEXT  NE 
NEXT  N 
REM  PRINT"  " 

REM  PRINT"  INPUT  OF  EXTERNAL  ELEMENTS" 

REM  PRINT"  " 

FOR  N-l  TO  NUMOS 

INPUT# 1,EXTRNL(N,1) ,EXTRNL(N,2) , 
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EXTRNL(N,3) , EXTRNL ( N , 4 ) , NUMMTX( N ) 

710  NEXT  N 

720  PRINT  #2,TITLE$ 

730  PRINT  #2,  "NUMBER  OF  STRAIGHT  LINE  SEGMENTS  TO  DEFINE 

BOUNDARY  ”,  NUMBS 

740  PRINT  #2, "NUMBER  OF  NON  BOUNDARY  POINTS  TO  CALCULATE 

RESULTS  AT  ",NUMOS 
750  IF  KSYM-1  THEN  GOTO  780 
760  IF  KSYM-2  THEN  GOTO  800 

770  IF  KSYM-3  THEN  GOTO  820  ELSE  GOTO  840 

780  PRINT  #2,"  "  \  PRINT  #2, "NO  SYMMETRY  CONDITIONS 
IMPOSED" 

790  GOTO  850 

800  PRINT  #2,"  "  \  PRINT  #2, "THE  LINE  X  -  XS  -  ";XSYM;"  IS 
A  LINE  OF  SYMMETRY" 

810  GOTO  850 

820  PRINT  #2,"  "  \  PRINT  #2, "THE  LINE  Y  -  YS  -  ";YSYM;"  IS 
A  LINE  OF  SYMMETRY" 

830  GOTO  850 

840  PRINT  #2,"  "  \  PRINT  #2, "THE  LINES  X  -  XS  -  ";XSYM;" 

AND  Y  -  YS  -  ";YSYM;"  ARE  LINES  OF  SYMMETRY" 

850  REM 

860  PRINT  #2,"  "  \  PRINT  #2 , "POISSON' S  RATIO  -  ";PR 
870  PRINT  #2,"  YOUNG'S  MODULUS  -  " ; E 

880  PRINT  #2,"  "  \  PRINT  #2 , "XX-COMPONENT  OF  FIELD  STRESS  - 
"  ;PXX 

890  PRINT  #2,"YY-COMPONENT  OF  FILED  STRESS  -  ";PYY 
900  PRINT  #2, "XY-COMPONENT  OF  FIELD  STRESS  -  ";PXY 
910  PRINT  #2,"  " 

920  PRINT  #2, "BOUNDARY  ELEMENT  DATA"  \  PRINT  #2,"  " 

930  PRINT  #2, "ELEMENT", "KODE","X  CENTER", "Y  CENTER" 

940  FOR  I-l  TO  NUMBE 

950  PRINT  #2,1, KOD(I) ,XM(I) ,YM(I) 

960  NEXT  I 

970  PRINT  #2,"  " 

980  PRINT  #2, "ELEMENT" , "LENGTH", "ANGLE" , "US  OR  SIGMA-S","UN 
OR  SIGMA-N" 

990  FOR  M-1  TO  NUMBE 
1000  MSIZE-2.0*A(M) 

1005  IF  COSBET(M)-0.0  AND  SINBET( M) >0 . 0  THEN  ANGLE-90  \  GOTO 
1020 

1007  IF  COSBET(M)-0.0  AND  SINBET(M)<0.0  THEN  ANGLE-270  \ 

GOTO  1020 

1010  ANGLE-180  *  ATN(SINBET(M)/COSBET(M) )/PI 

1015  IF  ANGLE<0  THEN  ANGLE-ANGLE+1 80 

1020  PRINT  #2,M,MSIZE,ANGLE,B( 2*M-1 ) ,B( 2*M) 

1030  NEXT  M 
1040  REM  PRINT"  " 

1050  REM  PRINT"  ADJUST  STRESS  BOUNDARY  VALUES  TO  ACCOUNT 
FOR  INITIAL  STRESSES  " 

1060  REM  PRINT"  " 

1070  FOR  N-1  TO  NUMBE 
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1080  NN-2*N 

1090  NS-NN-1 

1100  COSB-COSBET(N) 

1110  SINB»S1NBET(N) 

1120  SIGS-( PYY-PXX) *SINB*COSB+PXY* ( COSB*COSB-SINB*SINB ) 

1130  SIGN-PXX*SINB*SINB-2.0*PXY*SINB*COSB+PYY*COSB*COSB 

1140  IF  KOD(N)-l  THEN  GOTO  1170 

1150  IF  KOD(N)-2  THEN  GOTO  1240 

1160  IF  KOD(N)-3  THEN  GOTO  1200  ELSE  GOTO  1230 

1170  B(NS)-B(NS)-SIGS 

1180  B(NN)-B(NN)-SIGN 

1190  GOTO  1240 

1200  REM 

1210  B(NN)-B(NN)-SIGN 

1220  GOTO  1240 

1230  B(NS)-B(NS)-SIGS 
1240  NEXT  N 
1250  REM  PRINT"  " 

1260  REM  PRINT"COMPUTE  INFLUENCE  COEFFICIENTS  AND  SET  UP 
SYSTEM  OF  ALGEBRAIC  EQUATIONS" 

1270  REM  PRINT"  " 

1280  REM 

1290  FOR  I-l  TO  NUMBE 

1295  REM  PRINT"  "  \  REM  PRINT"  FOR  ELEMENT  ";I 
1300  IN-2*I 

1310  IS-IN-1 

1320  XI-XM(I) 

1330  YI-YM(I) 

1340  COSBI-COSBET( I ) 

1350  SINBI-SINBETd ) 

1360  KODE-KOD(I) 

1370  REM 

1380  FOR  J-1  TO  NUMBE 

1390  JN-2*J 

1400  JS-JN-1 

1410  REM  CALL  INITL 

1415  GOSUB  10000 

1420  XJ-XM(J) 

1430  YJ-YM(J) 

1440  COSBJ-COSBET( J) 

1450  SINBJ-SINBET( J) 

1460  AJ-A(J) 

1470  REM  CALL  COEFF ( XI , YI , XJ , YJ , AJ , COSB J , SINB J , +1 ) 

1480  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 

1490  QSIN-SINBJ  \  QQ-1 

1500  GOSUB  15000 

1510  IF  KSYM*1  THEN  GOTO  1690 

1520  IF  KSYM-2  THEN  GOTO  1550 

1530  IF  KSYM-3  THEN  GOTO  1580  ELSE  GOTO  1610 

1540  REM 

1550  XJ-2.0*XSYM-XM( J) 

1560  REM  CALL  COEFF ( XI , Y1 » XJ , YJ , AJ , COSB J , -SINB J , -1 ) 


1562  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
1564  QSIN— SINBJ  \  QQ— 1  \  GOSUB  15000 
1570  GOTO  1690 

1580  YJ-2.0*YSYM-YM( J) 

1590  REM  CALL  COEFF ( XI , YI , XJ , YJ ,AJ , -COSBJ , SINBJ , -1 ) 

1592  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS— COSBJ 
1594  QSIN-SINBJ  \  QQ— 1  \  GOSUB  15000 
1600  GOTO  1690 

1610  XJ-2.0*XSYM-XM( J) 

1620  REM  CALL  COEFF ( XI , YI , XJ , YJ , AJ , COSBJ , -SINBJ , -1 ) 

1622  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
1624  QSIN— SINBJ  \  QQ— 1  \  GOSUB  15000 
1630  XJ-XM(J) 

1640  YJ-2.0*YSYM-YM( J) 

1650  REM  CALL  COEFF ( XI , YI ,XJ , YJ ,AJ , -COSBJ , SINBJ , -1 ) 

1652  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \ 

QCOS-COSBJ 

1654  QSIN-SINBJ  \  QQ— 1  \  GOSUB  15000 
1660  XJ-2.0*XSYM-XM( J) 

1670  REM  CALL  COEFF ( XI , YI , XJ , YJ ,AJ , -COSBJ , -SINBJ , +1 ) 

1672  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \ 

QCOS-COSBJ 

1674  QSIN-SINBJ  \  QQ-1  \  GOSUB  15000 
1680  REM 
1690  REM 

1700  IF  KODE-1  THEN  GOTO  1740 

1710  IF  KODE-2  THEN  GOTO  1800 

1720  IF  KODE-3  THEN  GOTO  1860  ELSE  GOTO  1920 

1730  REM 

1740  C( IS, JS)-(SYYS-SXXS) *SINBI*COSBl 

+SXYS*(COSBI*COSBI-SINBI*SINBI ) 

1750  C( IS, JN)-( SYYN-SXXN)*SINBI*COSBl 

+SXYN*(COSBI*COSBI-SINBI*SINBI ) 

1760  C( IN, JS)-SXXS*SINBI*SINBI 

-2 . 0*SXYS*SINBI*COSBI+SYYS*COSBI*COSBI 
1770  C( IN, JN)-SXXN*SINBI*SINBI 

-2.0*SXYN*SINBI*COSBI+SYYN*COSBI*COSBI 
1780  GOTO  1970 

1790  REM 

1800  C( IS, JS)-UXS*COSBI+UYS*SINBI 
1810  C( IS, JN)-UXN*COSBI+UYN*SINBI 

1820  C( IN, JS)— UXS*SINBI+UYS*COSBI 

1830  C( IN, JN)— UXN*SINBI+UYN*COSBI 

1840  GOTO  1970 

1850  REM 

1860  C( IS, JS)-UXS*COSBI+UYS*SINBI 
1870  C( IS, JN)-UXN*COSBI+UYN*SINBI 

1880  C( IN, JS)-SXXS*SINBI*SINBI 

-2 .0*SXYS*SINBI*COSBI+SYYS*COSBI*COSBI 
1890  C( IN, JN)-SXXN*SINBI*SINB1 

-2.0*SXYN*SINBI*COSBI+SYYN*COSBI*COSBI 
1900  GOTO  1970 
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1910 

1920 

1930 

1940 

1950 

1970 

1975 

1980 

1990 

2000 

2010 

2020 

2030 

2040 

2050 

2060 

2069 

2070 

2080 

2090 

2092 

2100 

2110 

2120 

2130 

2140 

2150 

2160 

2170 

2180 

2190 

2200 

2210 

2220 

2230 

2240 

2250 

2260 

2270 

2280 

2290 

2300 

2310 

2320 

2330 

2340 

2350 


REN 

C(IS, JS)-( SYYS-SXXS)*SINBI*COSBI 
+SXYS*(COSBI*COSBI-SINBI*SINBI ) 

C ( IS , JN ) - ( SYYN-SXXN ) *SINBI *COSBI 
+SXYN*(COSBI*COSBI-SINBI*SINBI ) 

C( IN, JS)— UXS*SINBI+UYS*COSBI 
C( IN, JN)— UXN*SINBI+UYN*COSBI 
NEXT  J 
NEXT  I 

REM  PRINT"  " 

REM  PRINT"  SOLVE  SYSTEM  OF  ALGEBRAIC  EQUATIONS  " 
REM  PRINT"  " 

N-2*NUMBE 

REM  CALL  SOLVE(N) 

GOSUB  20000 
REM  PRINT"  " 

REM  PRINT"  COMPUTE  BOUNDARY  DISPLACEMENTS  AND 
STRESSES  " 

REM  PRINT"  " 

PRINT  #2,"  " 

PRINT  #2,"  DISPLACEMENTS  AND  STRESSES  AT  BOUNDARY 
ELEMENT  MIDPOINTS" 

PRINT  #2,"  " 

PRINT  #2 , "ELEMENT" , "UX" , "UY" , "US" , "UN" 

PRINT  #2,"SIGXX  SIGYY  SIGXY  SIGS 

SIGN  SIGT" 

FOR  I-l  TO  NUMBE 
XI-XM( I ) 

YI-YM( I ) 

COSBI-COSBET( I ) 

SINBI-SINBET( I ) 

REM 

UX-0.0 

UY-0.0 

SIGXX-PXX 

SIGYY-PYY 

SIGXY-PXY 

REM 

FOR  J-1  TO  NUMBE 
JN-2*J 
JS-JN-1 

REM  CALL  INITL 

GOSUB  10000 
XJ-XM( J) 

YJ-YM( J) 

AJ-A( J) 

COSBJ-COSBET{ J) 

SINBJ-SINBET( J) 

REM  CALL  COEFF(XI,YI,XJ,YJ,AJ,COSBJ,SINBJ,+l) 

QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
QSIN-SINBJ  \  QQ-1  \  GOSUB  15000 
IF  KSYM-1  THEN  GOTO  2650 
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2360  IF  KSYM-2  THEN  GOTO  2390 

2370  IF  KSYM-3  THEN  GOTO  2450  ELSE  GOTO  2510 

2380  REM 

2390  XJ-2.0*XSYM-XM( J) 

2400  REM  CALL  COEFF ( XI , YI , XJ , YJ , AJ , COSBJ , -SINE J , -1 ) 

2410  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 

2420  QSIN— SINBJ  \  QQ— 1  \  GOSOB  15000 

2430  GOTO  2650 

2440  REM 

2450  YJ-2.0*YSYM-YM( J) 

2460  REM  CALL  COEFF ( XI , YI , XJ , YJ ,AJ , -COSBJ , SINBJ , -1 ) 

2470  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \ 

QCOS— COSBJ 

2480  QSIN-SINBJ  \  QQ— 1  \  QQ— 1  \  GOSUB  15000 
2490  GOTO  2650 

2500  REM 

2510  XJ-2.0*XSYM-XM( J) 

2520  REM  CALL  COEFF ( XI , YI ,XJ , YJ ,AJ , COSBJ , -SINBJ , -1 ) 

2530  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
2540  QSIN— SINBJ  \  QQ— 1  \  GOSUB  15000 
2550  XJ-XM(J) 

2560  YJ-2.0*YSYM-YM( J) 

2570  REM  CALL  COEFF( XI , YI ,XJ , YJ ,AJ , -COSBJ , SINBJ , -1 ) 

2580  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS— COSBJ 
2590  QSIN-SINBJ  \  QQ— 1  \  GOSUB  15000 
2600  XJ-2.0*XSYM-XM( J) 

2610  REM  CALL  COEFF( XI , YI ,XJ , YJ ,AJ , -COSBJ , -SINBJ , +1 ) 

2620  QXI-XI  \  QYI-YI  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \ 

QCOS— COSBJ 

2630  QSIN-SINBJ  \  QQ-1  \  GOSUB  15000 
2640  REM 

2650  REM 
2660  REM 

2670  UX-UX+UXS*P( JS)+UXN*P( JN) 

2680  UY-UY+UYS*P( JS)+UYN*P( JN) 

2690  SIGXX-SIGXX+SXXS*P( JS)+SXXN*P( JN) 

2700  SIGYY-SIGYY-*-SYYS*P(  JS)+SYYN*P(  JN) 

2710  SIGXY-SIGXY+SXYS*P( JS)+SXYN*P( JN) 

2720  REM 

2730  NEXT  J 
2740  REM 

2750  US-UX*COSBI+UY*SINBI 

2760  UN— 1.0*UX*SINBI+UY*COSBI 

2770  SIGS-(SIGYY-SIGXX)*SINBI*COSBI 

+SIGXY*(COSBI*COSBI-SINBI*SINBI ) 

2780  SIGN-SIGXX*SINBI*SINBI 

-2 .0*SIGXY*SINBI*COSBI+SIGYY*COSBI*COSBI 
2790  SIGT-SIGXX*COSBI*COSBI 

+2 . 0*SIGXY*SINBI*COSBI+SIGYY*SINBI*SINBI 
2800  REM 

2810  OUTPT( I,1)-UX 
2820  OUTPT( I ,2)-UY 
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2830 

2840 

2850 

2860 

2870 

2880 

2890 

2900 

2905 

2910 

2912 

2920 

2930 

2935 

2936 
2940 
2950 

2960 

2965 

2990 

3000 

3010 

3020 

3030 

3040 

3042 

3045 

3050 

3060 

3070 

3080 

3090 

3100 

3110 

3120 

3130 

3140 

3150 

3160 

3170 

3180 

3190 

3200 

3210 

3220 

3230 

3240 


OUTPT( I,3)-US 
OUTPTd  ,4)-UN 
OUTPT(I,5)-SIGXX 
OUTPTd,  6  )-SIGYy 
OUTPT( I , 7 )-SIGXY 
OUTPTd, 8)-SIGS 
OUTPT( I ,9 )-SIGN 
OUTPT( I , 10 )-SIGT 

REM  PRINT"OUTPUT  FOR  ELEMENT  "  ;  I ;  "  COMPLETE" 

NEXT  I 

-  "  \  A$-A$+A$+A$+A$+A$+A$ 

FOR  I-l  TO  NUMBE 

PRINT  #2 , I ,OUTPT( 1,1) ,OUTPT( 1,2) ,OUTPT( I , 3 ) ,OUTPT( 1,4) 
PRINT  #2  USING  A$;  OUTPTd,  5),  OUTPTd,  6),  OUTPT  (1,7), 
OUTPTd,  8)  ,  OUTPTd,  9)  , OUTPTd,  10) 

PRINT  #2,  "  " 

NEXT  I 

REM  THIS  IS  THE  BEM  ELEMENT  DISP-STRESS  MATRIX  OUTPT  TO 
FILE 

MAT  PRINT  #3  ,  OUTPT 

MAT  PRINT  #3  ,  XM  \  CLOSE  #3 

REM  PRINT"  " 

REM  PRINT"  COMPUTE  DISPLACEMENTS  AND  STRESSES  AT 
SPECIFIED  POINTS  IN  THE  BODY" 

REM  PRINT"  " 

IF  NUMOS  <-  0  THEN  GOTO  3910 
PRINT  #2,"  "  \  PRINT  #2,"  " 

PRINT  #2,"  DISPLACEMENTS  AND  STRESSES  AT  SPECIFIED 
POINTS  IN  THE  BODY" 

PRINT  #2,"  "  \  PRINT  #2, "POINT" , "X  COORD" , "Y 
COORD" , "UX" , "UY" 

PRINT  #2,"  ", "SIGXX","SIGYY", "SIGXY"  \  PRINT  #2,"  " 

NPOINT-0 

FOR  N-1  TO  NUMOS 

XBEG-EXTRNL(N, 1 ) 

YBEG-EXTRNL(N,2) 

XEND-EXTRNL(N, 3 ) 

YEND-EXTRNL(N,4) 

NUMPB-NUMMTX(N) 

NUMP-NUMPB+1 
DELX- ( XEND-XBEG ) /NUMP 
DELY- ( YEND- YBEG ) /NUMP 
IF  NUMPB  >  0  THEN  NUMP-NUMP+1 
IF  (DELX"2+DELY''2)  -  0  THEN  NUMP-1 
REM 

FOR  NI-1  TO  NUMP 

XP-XBEG+ ( NI-1 ) *DELX 
YP-YBEG+ ( NI-1 ) *DELY 
REM 

UX-0.0 

UY-0.0 

SIGXX-PXX 
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3250 

3260 

3270 

3280 

3290 

3300 

3310 

3320 

3330 

3340 

3350 

3360 

3370 

3880 

3380 

3390 

3400 

3410 

3420 

3430 

3440 

3450 

3460 

3470 

3480 

3490 

3500 

3510 

3520 

3530 

3540 

3550 

3560 

3570 

3580 

3590 

3600 

3610 

3620 

3630 

3640 

3650 

3660 

3670 

3680 

3690 

3700 

3710 

3720 

3730 

3740 

3750 


SIGYY-PYY 

SIGXY-PXY 

REM 

FOR  J-1  TO  NUMBE 
JN-2*J 
JS-JN-1 

REM  CALL  INITL 

GOSUP  10000 

XJ-XM( J) 

YJ-YM( J) 

AJ-A( J) 

REM 

IF  SQR( (XP-XJ) "2+( YP-YJ) "2 )  <  (2.0*AJ)  THEN  GOTO 
REM 

COSBJ-COSBET( J) 

SINBJ-SINBET( J) 

REM  CALL  COEFF(XP, YP,XJ, YJ,AJ,COSBJ,SINBJ,+l ) 

QXl-XP  \  QYI-YP  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
QSIN-SINBJ  \  QQ-1  \  GOSUB  15000 
REM  GOTO  (840,810,820,830) ,KSYM 

IF  KSYM-1  THEN  GOTO  3750 
IF  KSYM-2  THEN  GOTO  3490 

IF  KSYM-3  THEN  GOTO  3550  ELSE  GOTO  3610 
REM 

XJ-2.0*XSYM-XM( J) 

REM  CALL  COEFF(XP, YP,XJ, YJ,AJ,COSBJ,-SINBJ,-l ) 

QXI-XP  \  QYI-YP  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
QSIN— SINBJ  \  QQ— 1  \  GOSUB  15000 
GOTO  3750 
REN 

YJ-2.0*XSYM-YM( J) 

REM  CALL  COEFF(XP, YP,XJ,YJ,AJ,-COSBJ, SINBJ, -1 ) 

QXI-XP  \  QYI-YP  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
QSIN-SINBJ  \  QQ— 1  \  GOSUB  15000 
GOTO  3750 
REM 

XJ-2.0*XSYM-XM( J) 

REM  CALL  COEFF(XP,YP,XJ,YJ,AJ,COSBJ, -SINBJ, -1 ) 

QXI-XP  \  QYI-YP  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
QSIN-SINBJ  \  QQ— 1  \  GOSUB  15000 
XJ-XM( J) 

YJ-2.0*YSYM-YM( J) 

REM  CALL  COEFF(XP, YP,XJ,YJ,AJ,-COSBJ, SINBJ, -1 ) 

QXI-XP  \  QYI-YP  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
QSIN-SINBJ  \  QQ— 1  \  GOSUB  15000 
XJ-2 . 0*XSYM-XM( J ) 

REM  CALL  COEFF(XP, YP,XJ,YJ,AJ,-COSBJ, -SINBJ, +1 ) 

QXI-XP  \  QYI-YP  \  QXJ-XJ  \  QYJ-YJ  \  QAJ-AJ  \  QCOS-COSBJ 
QSIN-SINBJ  \  QQ-1  \  GOSUB  15000 
REM 
REM 
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3760  REM 

3770  UX-UX+UXS*P( JS)+UXN*P( JN) 

3780  UY-UY+UyS*P( JS)+UYN*P( JN) 

3790  SIGXX-SIGXX+SXXS*P( JS)+SXXN*P( JN) 

3800  .  SIGYY-SIGYY+SYYS*P( JS)+SYYN*P( JN) 

3810  SIGXY-SIGXY+SXYS*P( JS)+SXYN*P( JN) 

3820  REM 

3830  NEXT  J 

3840  REN 

3850  NPOINT-NPOINT+1 

3860  PRINT  #2,NPOINT,XP,YP,UX,UY  \  PRINT  #2," 

” ,SIGXX,SIGYY,SIGXY 
3870  REM 

3880  NEXT  NI 

3890  NEXT  N 

3900  REM 

3910  REM 
3920  REM 

4000  GOTO  25000 

10000  REM  SUBROUTINE  INITL 

10010  REM 

10020  SXXS-0.0 

10030  SXXN-0.0 

10040  SYYS-0.0 

10050  SYYN-0.0 

10060  SXYS-0.0 

10070  SXYS-0.0 

10080  SXYN-0.0 

10090  REM 

10100  UXS-0.0 

10110  UXN-0.0 

10120  UYS-0.0 

10130  UYN-0.0 

10140  REN 

10150  RETURN 

15000  REM  SUBROUTINE  COEFF( X , Y, CX . CY , A , COSB , SINB , MSYM ) 

15010  REM 

15020  X-QXI  \  Y-QYI  \  CX-QXJ  \  CY-QYJ  \  A-QAJ  \  COSB-QCOS 

15030  SINB-QSIN  \  MSYM-QQ 
15040  REM 

15050  COS2B-COSB*COSB-SINB*SINB 

15060  SIN2B-2.0*SINB*COSB 

15070  REM 

15080  XB-(X-CX)*COSB+( Y-CY)*SINB 

15090  YB— 1 .0*(X-CX)*SINB+( Y-CY)*COSB 

15100  REM 

15110  R1S-(XB-A)*(XB-A)+YB*YB 

15120  R2S-(XB+A)*(XB+A)+YB*YB 

15130  FLl-.5*LOG(RlS) 

15140  FL2-.5*L0G(R2S) 

15150  FB2-CNST*(FL1-FL2) 

15160  IF  YB  <>  0  GOTO  15200 
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15170 

15180 

15190 

15200 

15210 

15220 

15230 

15240 

15250 

15260 

15270 

15280 

15290 

15300 

15310 

15320 

15330 

15340 

15350 

15360 

15370 

15380 

15390 

15400 

15410 

15420 

15430 

15440 

15450 

15460 

15470 

15480 

15490 

20000 

20010 

20020 

20030 

20040 

20050 

20060 

20070 

20080 

20090 

20100 

20110 

20120 


FB3-0 

IP  ABS(XB)  <A  THEN  FB3-CNST*PI 
GOTO  15210 

FB3--CNST*(ATN( ( XB+A)/YB ) -ATN( (XB-A)/YB) ) 
FB1-YB*FB3+CNST*( ( XB-A) *FLl- ( XB+A) *FL2 ) 

FB4-CNST* ( YB/R1S-YB/R2S ) 

FB5-CNST* ( ( XB-A ) /RlS- ( XB+A ) /R2S ) 

REM 

UXPS-COND* ( PR3*COSB*FBl+YB* ( SINB*FB2+C0SB*FB3 ) ) 
UXPN-COND* ( -PR3*SINB*FBl-YB* ( COSB*FB2-SINB*FB3 ) ) 
UYPS-COND* ( PR3*SINB*FB1-YB* ( C0SB*FB2-SINB*FB3 ) ) 
UYPN-COND* ( PR3*COSB*FBl-YB* ( SINB*FB2+COSB*FB3 ) ) 
REM 

SXXPS-FB2+PR2* ( C0S2B*FB2-SIN2B* FB3 ) 

+YB* ( C0S2B*FB4+SIN2B*FB5) 

SXXPN-FB3-PR1* ( SIN2B*FB2+C0S2B*FB3 ) 

+YB* ( SIN2B*FB4-C0S2B*FB5 ) 

SYYPS-FB2-PR2* ( C0S2B*FB2-SIN2B*FB3 ) 
-YB*(C0S2B*FB4+SIN2B*FB5) 

SYYPN-FB3+PR1* ( SIN2B*FB2+C0S2B*FB3 ) 

-YB*( SIN2B*FB4-C0S2B*FB5) 

SXYPS-PR2* ( SIN2B*FB2+C0S2B*FB3 ) 
+YB*(SIN2B*FB4-C0S2B*FB5) 

SXYPN-PRl* ( C0S2B*FB2-SIN2B*FB3 ) 
-YB*(C0S2B*FB4+SIN2B*FB5) 

REM 

UXS-UXS+MSYM*UXPS 

UXN-UXN+UXPN 

UYS-UYS+MSYM*UYPS 

UYN-UYN+UYPN 

REM 

SXXS-SXXS+MSYM*SXXPS 

SXXN-SXXN+SXXPN 

SYYS-SyYS+MSYM*SYYPS 

SYYN-SYYN+SYYPN 

SXYS-SXYS+MSYM*SXYPS 

SXYN-SXYN+SXYPN 

REM 

RETURN 

REM  SUBROUTINE  SOLVE(N) 

REM 

NB-N-1 

FOR  J-1  TO  NB 
L-J  +  1 

FOR  JJ=L  TO  N 

XM-C( JJ, J)/C( J, J) 

FOR  I-J  TO  N 

C( JJ,I)-C( JJ,I)-C( J,I)*XM 
NEXT  I 

B( JJ)-B( JJ)-B( J)*XM 
NEXT  JJ 
NEXT  J 
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20130 

20140 

20150 

20160 

20170 

20180 

20190 

20200 

20210 

20220 

20230 

20240 

25000 

25300 

25400 

25401 


REM 

P(N)-B(N)/C(N,N) 

FOR  J-1  TO  NB 
JJ-N-J 
L-JJ+1 
SUM-0 . 0 
FOR  I-L  TO  N 

SUM-SUM+C ( J J , I ) *P ( I ) 

NEXT  I 

P  {  J J ) - ( B ( J J ) -SUM ) /C ( J J , J J  > 
NEXT  J 
RETURN 

REM  PRINT" END  OF  PROCESSING" 
CLOSE  #1 
CLOSE  #2 
END 
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Appendix  C:  Computer  Program  TWOFS99_EX 

This  appendix  contains  the  listing  for  the  program 
TWOFS99_EX.  This  program  used  the  stress  versus  x  location 
output  from  TWOFS99  and  computed  the  stress  intensity  factor 
asa  function  of  distance,  r,  from  the  crack  tip.  The  equation 
used  was 

Rj  -  loyy  ,43) 

The  distribution  of  vs  r  was  only  taken  as  valid  from  a 
distance  five  to  ten  percent  of  the  crack  length  away  from 

the  crack  tip.  The  Kj  data  was  then  fit  through  linear 

2  2 

regression  analysis  against  r  .  The  rational  for  selecting  r 

over  an  r  distribution  is  explained  in  the  main  body  of  the 

text.  The  program  inputs  the  name  of  the  source  file,  the 

crack  length,  the  hole  diameter,  and  hole  pitch.  The  data 

that  fit  in  the  acceptable  distances  from  the  crack  tip  are 

printed  with  calculated  values,  and  the  final  regression 

fit  for  Kj  at  r-0  is  printed.  All  Kj  values  used  to  create 

the  0  factors  in  the  parametric  tension  strip  study  were 

calulated  by  this  program. 

The  program  is  written  in  VAX  BASIC  3.1  and  run  on  a 
VAX  8800. 
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1  DIM  OUTPT( 300,10) ,  XM(300),X(300>,SIGYy(300),R(300),K(300) 
10  PRINT  "  PROGRAM  TWOFS99_EX"  \  PRINT  "  " 

20  REM  TO  EXTRACT  DATA  FROM  OUTPT  FILES 

30  INPUT  "ENTER  OUTPT  FILE  NAME  ROOT  ";QIN$ 

31  INPUT  "ENTER  PITCH:  "; PITCH 

35  INPUT  "ENTER  CRACK  LENGTH  A  :  " ;A 

36  INPUT  "ENTER  HOLE  DIA:  ";DIA 

37  PRINT  "INPUT  FILE  ROOT:  ";QIN$ 

38  PRINT  "CRACK  LENGTH  A  : " ;A  \  PRINT  "HOLE  DIAMETER  :";DIA 

39  PRINT  "PITCH-  PITCH 

40  QOUT$  -  QIN$  +  ".OUTPT_EX" 

42  QIN$  -  QIN$  +  ".OUTPT" 

50  OPEN  QIN$  FOR  INPUT  AS  #1 

55  OPEN  QOUT$  FOR  OUTPUT  AS  #2 

56  PRINT  #2,"  PROGRAM  TWOFS99_EX"  \  PRINT  #2,"  " 

57  PRINT  #2,"  INPUT  FILE  :  ";QIN$ 

58  PRINT  #2,"  OUTPUT  FILE  :  ";QOUT$ 

59  PRINT  #2,"  "  \  PRINT  #2,"  CRACK  LENGTH  -  " ;A 

60  PRINT  #2, "HOLE  DIAMETER  :";DIA  \  PRINT  #2, "PITCH  -  ", 'PITCH 

\PRINT  #2,"  " 

61  FOR  I-l  TO  300 

62  FOR  J-1  TO  10 

64  INPUT  #1,  OUTPT(I,J) 

66  NEXT  J 
68  NEXT  I 
70  FOR  I-l  TO  300 
72  INPUT  #1,  XM(I) 

74  NEXT  I 

76  XMIN  -  DIA/2  +  A  +  0.05*A  \  XMAX  -  XMIN  +  0.05*A 
78  REM  CHECK  FOR  LONG  CRACK  PROBLEM 

80  IF  XMAX  <  PITCH-(DIA/2 )  THEN  GOTO  83 

81  XMAX  -  PITCH-( DIA/2)  \  XMIN  -XMAX  -  (  XMAX  -  DIA/2  - 

A)/2 

82  PRINT# 2, "LARGE  CRACK  WARNING" 

83  PRINT  #2, "XMIN  (5%  A)  -  ";XMIN  \  PRINT  #2, "XMAX  (10%  A)  - 

" ;XMAX 

85  KOUNT-0.0 

88  PRINT  #2,"  "  \  PRINT  #2,"  " 

89  PRINT  #2,  "ELEMENT",  "X  DIM",  "TIP  RAD",  "TIP  RAD  '“2",  "SIGMA 

YY" , "KI" 

90  FOR  I  -  2  TO  100 

91  REM  CHECK  BEM  2  TO  100 

92  REM  CHECK  FOR  5%  <  X  <  10%  OF  A 

94  IF  XM(I)  >  XMAX  OR  XM{I)  <  XMIN  THEN  GOTO  180 
100  KOUNT-KOUNT  +  1 

110  X(I)  -  XM(I)  \  SIGYY(I)  -  OUTPT( I ,6)/1000 
120  R(I)  -  X(I)  -  A  -  DIA/2. 0 
125  R2-  R(I)'2 

130  K(I)  -  SIGYY(I)  *  (  2  *  PI  *  R(I)  )''0.5 
140  SUMR  -  SUMR  +  R2 
150  SUMR2  -  SUMR2  +  R2'2 
160  SUMRK  -  SUMRK  +  R2*K(I) 
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170  SUMK  -  SUMK  +  K(I) 

175  PRINT  #2,  I,  XM(1),  R(I),R(I)^2,  SIGYY(I),  K(I) 

180  NEXT  I 

200  B  -  (KOUNT  *  SUMRK  -  SUMR*SUMK )/( KOUNT  *SUMR2- { SUMR) " 2 ) 

210  KICFIT  -  (SUMK  -  B*SUMR )/KOUNT 

215  BETA  -  KICPIT  /  (  46  *  SQR(  PI  *  A)  ) 

220  PRINT  #2  - " 

230  PRINT  #2  " 

240  PRINT  #2  KI  (REGRESSION  FIT  R-0.0)  -  ";KICFIT 
245  PRINT  #2  BASED  ON  R  SQUARED  " 

250  PRINT  #2, "  " 

255  PRINT  #2,"  BETA  (SlG-46)  -  ";BETA 

260  PRINT  #2,  '• - " 


4 
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Appendix  D:  Computer  Program  CHOLE 


This  appendix  contains  the  listing  for  the  program 
CHOLE.  This  program  is  a  model  generator  for  the  tension 
strip  parametric  study  of  section  VII.  The  input  to  the 
program  is  hole  diameter,  pitch,  and  crack  length.  The 
program  divides  the  crack  into  segments  with  the  F.R.  Harris 
refinement  technique  [6].  The  final  model  as  output  is  in  a 
format  required  for  TWOFS99  to  read  in. 

All  of  the  models  used  in  the  tension  strip 
parametric  study  were  created  with  CHOLE.  CHOLE  is  a  VAX 
BASIC  3.1  program  run  on  a  VAX  8800. 


20  PRINT  "BEM  HOLE  WITH  CRACK  MODEL  GENERATOR  -  QUAD  FINITE 
BOUND" 

22  print  "  GRADUATED  CRACK  ELEMENTS  3-3-3-25  RULE  " 

30  PRINT  "  " 

35  INPUT"ENTER  NAME  OF  OUTPUT  FILE  :  ";0$ 

36  OPEN  0$  FOR  OUTPUT  AS  #1 

40  INPUT"ENTER  HOLE  DIAMETER" ;DIA 

50  INPUT"ENTER  DISTANCE  BETWEEN  HOLE  CENTERS  ", •PITCH 
60  INPUT"ENTER  LENGTH  OF  CRACK  " ;A 

70  PRINT  #1,"  TWO  HOLES  S-46  D-";DIA;"  P-" ; PITCH ; "  A-";A 

80  SXX-0.0  \  SYY-46000.  \  SXY-0.0 

115  RAD  -  DIA/2.0 

120  CIRCUM-  2  *  3.14159  *  RAD 

130  REM  DIVIDE  CRACK  BY  20  TO  GET  ELEMENT  LENGTH 
140  ELEN  -  A/12 

150  REM  CALCULATE  HOW  MANY  ELEMENTS  IN  HALF  CIRCLE  (HOLE) 

160  CEL  -(  CIRCUM/  ELEN  )/2.0 

170  CEL  -  INT(  CEL)  +1  \  IF  CEL  <  20  THEN  CEL-20 

175  REM  4  FOR  CRACK  4  FOR  PRECRACK  1  FOR  INBETWEEN 

176  REM  2  CLOSE  HORIZON  2  SIDES  1  TOP  1  SPC 
185  ELTOT  -  2  *  CEL  +  4  +  44-1  +  2  +  2  +  1  +  1 
187  IF  RAD  +  A  -  PITCH/2  THEN  ELTOT-ELTOT-1 
190  PRINT  #1,  ELTOT 0,1, .3,10. 3E6" 

220  PRINT  #1,"0.0" 

230  PRINT  #1, "0.0" 

240  PRINT  #1,"0.0" 

242  T$  -  "1,0,0"  \  C$-"," 
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244  REM  THIS  IS  THE  NON  CRACK  MATERIAL  BETWEEN  HOLES 

245  IF  PITCH  >  2  *  (  A  +  RAD)  THEN  GOTO  270 

247  LTEMP  -  PITCH  -  PAD  -  RAD  -  A 

248  Xl-PITCH  -  RAD  \  X2  -  A  +  RAD  +  0.5  *  LTEMP 


249 

PRINT 

#1,  "15," 

• 

9 

Xl  ;  C$  ;  Yl  ;  C$  ; 

X2 

• 

9 

c$ 

9 

9 

Y2 

9 

9 

c$ 

• 

"4,0,0" 

250 

X1-X2 

\  X2  -  A 

+ 

RAD  +  0.25  *  LTEMP 

251 

PRINT 

#1,  "15," 

! 

Xl  ;  C$  ;  Yl  ;  C$  ; 

X2 

9 

9 

c$ 

• 

9 

Y2 

• 

9 

c$ 

? 

"4,0,0" 

252 

X1-X2 

\  X2  -  A 

+ 

RAD  +  0.125  *  LTEMP 

253 

PRINT 

#1,  "15," 

• 

9 

Xl  ;  C$  ;  Yl  ;  C$  ; 

X2 

• 

9 

c$ 

9 

9 

Y2 

9 

9 

c$ 

9 

"4,0,0" 

254 

X1-X2 

\  X2  -  A 

+ 

RAD 

255 

PRINT 

#1,  "15," 

« 

9 

Xl  ;  C$  ;  Yl  ;  C$  ; 

X2 

9 

9 

c$ 

• 

9 

Y2 

9 

9 

c$ 

• 

9 

"4,0,0" 

256 

X1-X2 

\  X2  -  A 

+ 

RAD  -  0.125  *  LTEMP 

257 

PRINT 

#1,  "25," 

« 

9 

Xl  ;  C$  ;  Yl  ;  C$  ; 

X2 

• 

9 

c$ 

9 

9 

Y2 

9 

9 

c$  ; 

T$ 


258 

X1-X2  \  X2 

-  A  -t-  RAD 

-  0.25  * 

LTEMP 

259 

PRINT  #1, 

"3,";  Xl  ; 

C$  ;  Yl 

;  C$  ; 

X2 

;  c$ 

;  y2 

;  c$ 

• 

9 

T$ 

260 

X1-X2  \  X2 

-  A  +  RAD- 

0.5  * 

LTEMP 

261 

PRINT  #1, 

" 3 , " ;  Xl  ; 

C$  ;  Yl 

;  C$  ; 

X2 

;  c$ 

;  y2 

;  c$ 

9 

9 

T$ 

262 

X1-X2  \  X2 

-  A+  RAD  - 

•  LTEMP 

263 

PRINT  #1, 

"3,";  Xl  ; 

C$  ;  Yl 

;  c$  ; 

X2 

;  C$ 

;  Y2 

;  c$ 

9 

9 

T$ 

264  IF  A+RAD  -  PITCH/2  THEN  GOTO  290 

265  LTEMP  -  X2  -  RAD  \  LTEMP2  -  (XI  -  X2)/3 

266  LTOT  -  INT( LTEMP/LTEMP2 )  +  1 

267  X1-X2  \  X2  -  RAD  \  Y2  -  0.  \  Yl  -  0. 

268  PRINT  #l,LTOT  XI  ;  C$  ;  Yl  ;  C$  ;  X2  ;  C$  ;  Y2  ;  C$ 

;  "  1 , 0 , 0  " 

269  GOTO  290 

270  LTEMP  -  PITCH-RAD  -RAD  -A  -A 

271  LTOT  -  INT  (LTEMP/A)  +1 

272  XI  -  PITCH-RAD  \  Yl  -  0 . 0  \  X2  -  XI  -  LTEMP  \  Y2  -  Yl 

273  PRINT  #l,LTOT  XI  ;  C$  ;  Yl  ;  C$  ;  X2  ;  C$  ;  Y2  ;  C$ 

;"4,0,0" 

274  X1-X2  \  X2  -  A  +  RAD  +  0.5  *  A 

275  PRINT  #1,  "3,";  Xl  ;  C$  ;  Yl  ;  C$  ;  X2  ;  C$  ;  Y2  ;  C$ 

; " 4 , 0 , 0" 

276  X1-X2  \  X2  -  A  +  RAD  +  0.25  *  A 

277  PRINT  #1,  "3,";  Xl  ;  C$  ;  Yl  ;  C$  ;  X2  ;  C$  ;  Y2  ;  C$ 

;"4,0,0" 

278  X1-X2  \  X2  -  A  +  RAD  +  0.125  *  A 

279  PRINT  #1,  "3,”;  Xl  ;  C$  ;  Yl  ;  C$  ;  X2  ;  C$  ;  Y2  ;  C$ 

."4,0,0" 

280  X1-X2  \  X2  »  A  +  RAD 

281  PRINT  #1,  "25,";  Xl  ;  C$  ;  Yl  ;  C$  ;  X2  ;  C$  ;  Y2  ;  C$ 

; "4 , 0 , 0" 

282  X1-X2  \  X2  -  A  +  RAD  -  0.125  *  A 
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283  PRINT  #1,  ’’25,";  XI  ;  C$  ;  Yl  ;  C$  ;  X2  ;  C$  ;  Y2  ;  C$  ; 
T$ 


284 

X1-X2  \  X2 

- 

A  +  RAD 

-  0.25 

*  A 

285 

PRINT  #1, 

"3 

"  •  Xl  • 

C$  ;  Yl 

;  c$ 

;  X2 

;  c$ 

;  Y2 

;  c$ 

• 

i 

T$ 

286 

X1-X2  \  X2 

- 

A  +  RAD- 

0.5  * 

A 

287 

PRINT  #1, 

"3 

"  •  Xl 

C$  ;  Yl 

;  c$ 

;  X2 

;  c$ 

;  Y2 

;  c$ 

• 

T$ 

288 

X1-X2  \  X2 

- 

RAD 

289 

PRINT  #1, 

"3 

"  •  XI 
/  /  / 

C$  ;  Yl 

;  c$ 

;  X2 

;  c$ 

;  Y2 

;  c$ 

• 

9 

T$ 


290 

291 
310 
320 
330 
340 
350 
400 

420 

600 

605 

607 

610 

620 

630 

640 

650 

660 

720 

721 

722 
728 

732 

734 

738 

750 

752 

754 

756 

758 

759 

760 
765 
770 
775 


REM  THIS  IS  THE  HOLE  CALCULATION  SECTION 

ANGLE  -  3.14159  \  DELA  -  ANGLE/  CEL  \  ANGLE-0.0 

FOR  I-l  TO  CEL 

XI  -  X2  \  Yl  -  Y2 

ANGLE  -  ANGLE  +  DELA 

X2  -  RAD  *  COS(  ANGLE) 

Y2  »  RAD  *  SIN(  ANGLE) 

PRINT  #1,"1,"  ;  XI  ;  C$  ;  Yl  ;  C$  ;  X2  ;  C$  ;  Y2  ;  C$  ; 
T$ 

NEXT  I 

REM  THIS  IS  THE  SECOND  HOLE 

ANGLE  -  3.14159  \  DELA  -  ANGLE/  CEL  \  ANGLE-0.0 
X2  -  RAD  +  PITCH\  Y2  -  0.0  \  C$-"," 

FOR  I-l  TO  CEL 

XI  -  X2  \  Yl  -  Y2 

ANGLE  -  ANGLE  -t-  DELA 

X2  -  RAD  *  COS(  ANGLE)  +  PITCH 

Y2  -  RAD  *  SIN(  ANGLE) 

PRINT  #l,’’l,"  ;  XI  ;  C$  ;  Yl  ;  C$  ;  X2  ;  C$  ;  Y2  ;  C$  ; 
T$ 

NEXT  I 

REM  THIS  IS  3-D  ON  LEFT  OF  LEFT  HOLE 

DIST-3*DIA  \X1  -  -RAD  \  X2  -  Xl  -  DIST  \  Yl-0.0  \  Y2-0.0 
PRINT  #1,"10,’’  ;  XI  ;  C$ ;  Yl  ;  C$ ;  X2 ;  C$ ;  Y2 ;  C$  ; 
"4,0,0" 


REM  THIS  IS  3-D  ON  RIGHT  OF  RIGHT  HOLE 

Xl  -  PITCH  +  RAD  +DIST  \  X2  -  Xl  -  DIST  \  Yl-0.0  \  Y2-0.0 
PRINT  #1,"10,’’  ;  Xl  ;  C$ ;  Yl  ;  C$;  X2 ;  C$ ;  Y2 ;  C$; 

"4,0,0" 


REM  THIS  IS  THE  FINITE  ( 3-DIA)  BOUNDARY 
REN  L  SIDE 

Xl  -  -RAD  -  DIST  \  X2  -  Xl 
Y2  -  DIST  +  RAD  \  Yl  -  0.0 

PRINT  #1,"10,’’  ;  Xl  ;  C$;  Yl  ;  C$;  X2 ;  C$ ;  Y2;  C$ ;  "1, 
0  ,  ;  SXX 
REM  TOP 
DIST  -  3  *  DIA 

Xl  -  -RAD  -  DIST  \  X2  -  RAD  +  PITCH  +  DIST 
Yl  -  DIST  +  RAD  \  Y2  «  Yl 

PRINT  #1,"40,’’  ;  Xl  ;  C$;  Yl  ;  C$ ;  X2 ;  C$ ;  Y2 ;  C$ ;  "1, 
0  ,  ;  SYY 
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780  REM  R  SIDE 
782  XI  -  X2 

784  Yl  -  DIST  +  RAD  \  Y2  -  0.1 

785  PRINT  #1,"10,"  ;  Xl  ;  C$;  Yl  ;  C$;  X2 ;  C$ ;  Y2 ;  C$ ;  "1, 

0 , " ; SXX 

790  REM  THIS  IS  THE  SPC 
800  Y1-Y2  \  Y2  -  0.0 

815  PRINT  #l,''l,"  ;  Xl  ;  C$ ;  Yl  ;  C$;  X2 ;  C$;  Y2 ;  C$ ;  "2, 

0,0" 

900  PRINT  "END  OF  PROCESSING" 

910  CLOSE  #1 
1000  END 
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Appendix  E:  Fitting  NASA/FLAGRO  Crack  Growth  Output 


This  appendix  contains  the  output  from  the 
NASA/FI -'.GRO  analysis  of  section  VIII.  All  three  of  the 
analysis  used  the  same  materials  and  stress  spectrums. 


A.  Bowi;.-  Solution  Analysis 

FAilGUE  CRACK  GROWTH  ANALYSIS 

(computed:  NASA/FLAGRO,  1986  Aug  version,  1987  Jul  rev.) 
U.S.  customary  units  (inches,  ksi,  ksi  sqrt(in)l 

PROBLEM  TITLE 

TEST  OF  BOWIE  SOLUTION  ANALYSIS 
GEOMETRY 

MODEL:  TC03-Through  crack  from  hole  in  plate. 

Plate  Thickness,  t  -  0,2500 

"  Width,  W  -  100.0000 

Hole  Diameter,  D  -  0.2500 

Distance  of  Hole  Center  to  Edge,  B  50.0000 

FLAW  SIZE: 

a  (init.)  •  0.7500E-01 


MATERIAL 


MATL  1:  7075-T6  AL,  L-T 
Material  Properties: 


Matl 

No. 

YS  : 

Kle  : 

Klc  : 

Ak 

Bk 

Thk  :  Kc  : 

KIscc : 

1 

65.  o’: 

42.0*: 

27.0: 

0.75: 

1.25 

0.250:  54.9: 

Matl 

Crack  Growth  Eqn  Constants  (closure) 

No. 

C 

:  n 

:  p  : 

q  : 

DRo 

Co  :  d  :  DKl 

: Alpha : Smax/ 

: 

:  :SIGo 

1 

0.275D-07:2. 836:0. 50: 

0.50i 

2.50 

1.00:1.00:  5.74 

':  1.75  5  0.30 
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TEST  OF  BOWIE  SOLUTION  ANALYSIS 
MODEL:  TC03 

FATIGUE  SPECTRUM  STRESS  TABLE 


S  : 

M 

NUMBER 

so 

SI 

T  :  A 
E  ;  T 

OF 

FATIGUE 

( ksi  ) 

(  ksi  ) 

P  : 

L 

CYCLES 

tl  :  t2 

tl  :  t2 

1 : 

1 

100 

0.00:  30.00 

o 

1  o 

1  • 

1  o 

1 

1 

1 

1  o 

1  o 
1  • 

1  o 

1 

1 

Environmental  Crack  Growth  Check  for  Sustained  Stresses 
( Kmax  less  than  KIscc):  NOT  SET 


TEST  OF  BOWIE  SOLUTION  ANALYSIS 
MODEL:  TC03 

ANALYSIS  RESULTS: 


Block 

Final  Flaw  Size 

K  max 

Step 

a 

a-tlp 

1 

0.092832 

23.873066 

2 

0.111912 

24.298998 

3 

0.132318 

24.740814 

4 

0.154225 

25.221426 

5 

0.177883 

25.756659 

6 

0.203603 

26.360479 

7 

0.231775 

26.993177 

8 

0.262879 

27.740626 

9 

0.297515 

28.528124 

10 

0.336449 

29.429157 

11 

0.380675 

30.424570 

12 

0.431521 

31.576437 

13 

0.490817 

32.815969 

14 

0.561179 

34.270541 

15 

0.646548 

35.996892 

16 

0.753263 

38.005786 

17 

0.892600 

40.549542 

18 

1.088274 

43.766852 

19 

1.412186 

48.761760 

FINAL  RESULTS: 

Unstable  crack  growth,  max  stress  intensity  exceeds  critical  value 
K  max  ■  55.00  K  cr  -  54.94 

at  Cycle  No.  56.  of  Load  Step  No.  1  of  Block.  No.  20 

Crack  Size  a  -  1.87187 
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B .  Sh\vakumar  Solution  Analysis 


FATIGUE  CRACK  GROWTH  ANALYSIS 


(compu*’ed:  NASA/FLAGRO,  1986  Aug  version,  1987  Jul  rev.) 
U.S.  customary  units  [inches,  )tsi,  ksi  sqrt(in)) 

PROBLEM  TITLE 


TEST  OF  SHIVAKUMAR  SOLUTION  ANALYSIS 
GEOMETRY 


MODEL;  TC05-Through  crack  from  hole  in  row  of  holes. 

Plate  Thickness,  t  -  0.2500 

Hole  Diameter,  D  ••  0.2500 

Distance  between  Holes,  H  >  1.0000 

Ratio  of  Hole  Diameter  to  Edge  Distance,  D/B  -  0.0000 

(Ratio  of  0.0  denotes  a  very  large  edge  distance) 

FLAW  SIZE: 

a  (init.)  -  0.7500E-01 


MATERIAL 


MATL  1:  7075-T6  AL,  L-T 
Material  Properties: 


Matl 

YS  ; 

Kle 

Klc 

Ak 

:  Bk 

Thk 

Kc  : 

KI  see : 

No . 

1 

65.0: 

42  . 

0; 

27.0 

0.7 

5:  1.25 

0.2 

50: 

54.9: 

Matl 

Crack  Growth 

Eqn  Constants 

( closure ) 

No. 

C 

n 

P 

q 

;  DKo 

Co  : 

d 

DKl 

: Alpha 

Smax/ 

SIGO 

1 

0.275D- 

07:2. 

836 

0.50 

0.50 

i  2.50 

1.00: 

1.00 

5.74 

i  1.75 

0.30 
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TEST  CF  SHIVAKUMAR  SOLUTION  ANALYSIS 
MODEL;  TC’S 

FATIGUE  SPECTRUM  STRESS  TABLE 


SAWTOOTH  0  -  30  KSI 


S  :  M 

NUMBER 

SO 

SI 

S2 

T  :  A 

OF 

E  :  T 

FATIGUE 

(  ksi  ) 

( ksi  ) 

(ksi) 

P  :  L 

CYCLES 

tl  ;  t2 

tl  :  t2 

tl  :  t2 

1:  1 

100 

0.00:  30.00 

1 

1 

1 

O  1 
•  1 
O  1 
O  1 

1 

1 

1 

O  1 
•  1 
O  1 
O  1 

1 

1 

O  1 
.  1 
O  1 
O  1 

1 

1 

1 

O  1 
.  1 
O  ( 
O  1 

Environmental  Crack  Growth  Check  for  Sustained  Stresses 
( Kmax  less  than  KIscc);  NOT  SET 


TEST  OF  SHIVAKUMAR  SOLUTION  ANALYSIS 
MODEL:  TC05 

ANALYSIS  RESULTS: 


ADVISORY:  Estimated  Net  Section  Stress  > 

Yield  Strength. 

at  Cycle  No. 
Crack  Size 

0.  of  Load  Step  No. 
a  -  0.750000E-01 

1  of  Block  No. 

Block 

Final  Flaw  Size 

K  max 

Step 

a 

a-tip 

1 

0.094329 

24.414055 

2 

0.115210 

24.922510 

3 

0.137799 

25.448033 

4 

0.162406 

26.050727 

5 

0,189473 

26.739437 

6 

0.219583 

27.535714 

7 

0.253507 

28.430209 

8 

0.292277 

29.476453 

9 

0.337329 

30.668877 

10 

0.390770 

32.098262 

11 

0.456008 

33.810283 

12 

0.539788 

36.204720 

13 

0.668205 

42.025305 

FINAL  RESULTS:  ,  , 

Unstable  crack  growth,  max  stress  intensity  exceeds  critical  value 
K  max  -  56.32  K  cr  -  54.94 

at  Cycle  No.  24.  of  Load  Step  No.  1  of  Block  No.  14 

Crack  Size  a  «  0.747562 
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C .  g  Look-Up  Table  Analysis 


FATIGUE  CRACK  GROWTH  ANALYSIS 

(computed:  NASA/FLAGRO,  1986  Aug  version,  1987  Jul  rev.) 
U.S.  customary  units  [inches,  )(si,  ksi  sqrt(in)l 

PROBLEM  TITLE 


TEST  OF  BOUNDARY  ELEMENT  LOOK  UP  TABLE  ANALYSIS 
GEOMETRY 


MODEL:  DTOl-One-dimensional  data  table  for  through  crack. 
Plate  Thickness,  t  -  0.2500 

a/D  :  FO 


0.1000 
0.2000 
0.3000 
0.4000 
0.5000 
0.6000 
0.7000 
0.8000 
0.9000 

where 

SO  ;  TENSION  STRESS 

FLAW  SIZE; 

a  (init.)  -  0.7500E-01 


1.7400 
1 . 3800 
1.2500 
1 . 2000 
1.1900 
1.2000 
1.2400 
1.3400 
1.6800 


MATERIAL 


MATL  1:  7075-T6  AL,  L-T 
Material  Properties; 


Matl 

No. 

YS  : 

Kle  : 

Klc 

Ak 

Bk 

Thk  : 

Kc  : 

KIscc : 

1 

65. oi 

4  2.0: 

27.0 

0.75 

1.25 

0.250: 

54.9: 

Matl 

Crack  Growth  Eqn  Constants  (closure) 

No. 

C 

n 

:  P 

q  ’• 

DKo 

Co  :  d 

DKl 

: Alpha : Smax/ 

; 

:SIGo 

1 

0.275D-07:2.836;0.50 

0.50; 

2.50 

o 

o 

o 

o 

tH 

5.74 

:  1.75;  0.30 
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TEST  OF  BOUNDARY  ELEMENT  LOOK  UP  TABLE  ANALYSIS 
MODEL:  DTOl 


FATIGUE  SPECTRUM  INPUT  TABLE 


SAWTOOTH  0  -  30  KSI 

[Note:  Stress  -  Input  Val  .e  *  Stress  Factor] 
Stress  Factor  SFO:  1.00 


S  :  M 

NUMBER 

SO 

T  :  A 

OF 

E  :  T 

FATIGUE 

P  :  L 

CYCLES 

tl  :  t2 

1 :  1 

100 

0.00:  30.00 

Environmental  Crack  Growth  Check  for  Sustained  Str 
( Kmax  less  than  KIscc);  NOT  SET 


TEST  OF  BOUNDARY  ELEMENT  LOOK  UP  TABLE  ANALYSIS 
MODEL:  DTOl 

ANALYSIS  RESULTS: 


Block 

Step 

1 

2 

3 

4 

5 

6 


Final  Flaw  Size 
a 

0.100650 

0.131098 

0.166947 

0.211160 

0.270492 

0.361315 


K  max 

a-tip 

26.642060 

27,751991 

29.046889 

30.923850 

33.563043 

38.080432 


FINAL  RESULTS: 

Unstable  crack  growth,  max  stress  intensity  exceeds  critical  value 
K  max  -  55.11  K  cr  -  54.94 

at  Cycle  No.  95.  of  Load  Step  No.  1  of  Block  No.  7 

Crack  Size  a  -  0.599436 
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Appendix  F:  Comparison  of  Regression  Fit  Analysis 


This  appendix  shows  examples  of  linear  regression 

fits  for  large  and  small  crack  ratios,  for  both  r  (distance 

2 

from  crack  tip)  and  r  .  As  was  explained  in  the  main  text  of 

this  thesis,  there  is  insignificant  differences  for  the 

values  of  predicted  for  small  crack  ratio  problems  from 

2 

linear  regression  fits  of  vs  r  or  r  .  Figures  27  and  28 

2 

show  the  plots  of  vs  r  and  r  respectively  for  a  small 

crack  ratio  problem  (crack  ratio  -  0.1)  from  the  tension 

strip  parametric  study  of  section  VII.  This  data  is  for  hole 

diameter  equal  to  0.25  inches  and  pitch  equal  to  four 

diameters.  Both  plots  indicate  a  value  of  Kj  at  r-0  of 

1/2 

approximately  39  KSI(in)  .  However,  Figures  29  and  30  show 

the  same  plots  for  a  crack  ratio  of  0.9  (same  diameter  hole 

2 

and  pitch).  While  the  r  regression  fit  will  indicate  a  Kj 
value  of  112.7  KSI(in)^^^,  the  r  fit  data  will  not  even 
predict  a  positive  value  of  Kj.  It  is  hypothesized  that  the 
indicated  values  of  are  not  linear  in  r,  and  the  portion 
of  the  Kj  vs  r  curve  plotted  in  Figure  29  is  quadratic  in  r. 
Therefore  a  linear  regression  fit  is  inadequate  for  the  large 
crack  ratios,  and  was  not  used  for  the  parametric  study  of 
section  VII. 
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Figure  28.  Stress  Intensity  Factor  Vs  Radius"2  (Crack  Ratio-0.1) 
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Figure  29.  Stress  Intensity  Factor  Vs  Radius  (Crack  Ratio-0.9) 
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Figure  30.  Stress  Intensity  Factor  Vs  Radius''2  (Crack  Ratio“0i9) 
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